Understanding the dynamics of neuronal activation in response to sensory stimuli is pivotal for unraveling the complexities of information processing in the brain. This exploration extends beyond mere activation, focusing on identifying stimuli that elicit optimal neuronal responses, thereby maximizing their firing rates. \
While Linear-Nonlinear (LN) models have traditionally demonstrated a significant predictive power for responses in linear systems such as the retina, they show a lower predictive efficacy in identifying optimal sensory inputs for neurons characterized by nonlinear sensitivity. This complexity is particulary evident in primary visual cortex (V1), especially in response to naturalistic stimuli, as this involves navigating an intractably high-dimensional space of potential stimuli.
In response to these challenges, model-driven stimulus optimization leverages functional models capable of accurately predicting neuronal responses to a broad range of stimuli, including natural scenes. This promising approach facilitates an unrestricted exploration through a high-dimensional stimulus space, aiming to identify inputs that maximally activate specific neurons.
This research is committed to employing advanced deep learning techniques to identify and synthesize the Most Excitatory Input (MEI) for a single neuron in the visual system of a mouse, specifically through the use of end-to-end trained convolutional neural networks (CNNs). \
Most Excitatory Inputs (MEIs) represent the specific visual stimuli that are designed to maximally activate the neurons. These in silico stimuli are pivotal for delineating the functional characteristics of a neuronal single-cell, highlighting the features within a neuron's receptive field to which it is most responsive.
By adopting this innovative methodology, we aim to construct a stimulus that successfully combining the stimuli's features, is most likely to trigger the highest firing rates in neurons. Finally, we seek to investigate some properties of these artificially optimized stimuli and try to test the robustness of our findings using Gabor filters as benchmarks.
Our study leverages a highly specialized dataset provided by The Allen Institute Brain Observatory, which comprises neuronal response recordings from mice observing visual stimuli carefully adjusted for eye shape. The recording of these responses utilizes high-density electrodes stacked upon one another, capable of measuring field potential changes at a high rate. \ \
Here is an outline of the approach followed for this study:
Filtering and Selection of Neurons: The study begins with a meticulous process of filtering to select the most informative neurons from the recorded sessions. This selection is based on the neurons’ responsiveness and their potential to yield significant insights regarding the stimuli they process. This ensures the focus is maintained on data that is most likely to enhance our understanding. \
Training the CNN on Individual Neurons: For each of the selected neurons, a convolutional neural network is trained using the extensive dataset of natural scenes provided by The Allen Institute. This training is aimed at fine-tuning the neural network to predict the firing patterns of each neuron in response to diverse visual inputs, thereby optimizing the network’s ability to model individual neuronal responses. By training the CNN on these neural response data, the model not only predicts but will also aid in optimizing stimuli to significantly enhance the activation of specific neurons. \
Generation of Most Excitatory Inputs (MEIs): With the trained CNNs, the next step involves generating MEIs. This is done through an optimization process where initial stimuli are iteratively refined to increase the neuron's firing rate, under the guidance of the CNN’s predictive capabilities. \
Analysis and Comparison with Gabor Filters: After generating the MEIs, the study analyzes these inputs to identify the critical features and patterns that effectively stimulate neuronal activity. These results are then compared with responses induced by traditional Gabor filters, which are renowned for their efficacy in simulating responses of simple cells in the visual cortex. \
Conclusions and Implications: The study culminates in drawing conclusions from the comparative analysis between the CNN-generated MEIs and the Gabor filter responses. These insights aim to deepen our understanding of how neurons encode visual information and discuss limitations and insights for further analysis.
import pandas as pd
from scipy.stats import entropy
import numpy as np
# For visualisation
import seaborn as sns
import matplotlib.pyplot as plt
# For Anomaly Detection
from sklearn.ensemble import IsolationForest
# For Oracle Correlation
from scipy.stats import pearsonr
import random
# For Rasterplot
from rastermap import Rastermap
from scipy.stats import zscore
# For AllenSDK
import os
from allensdk.brain_observatory.ecephys.ecephys_project_cache import EcephysProjectCache
# Filter Ignore
# import warnings
# warnings.filterwarnings('ignore')
from collections import OrderedDict, defaultdict
import numpy as np
import torch
import torch.nn as nn
from itertools import product
from torch.nn import functional as F
import numpy as np
import scipy.signal
import torch
import math
from math import ceil
from torch import nn as nn
# from .module import Module
from torch.nn import Parameter
from torch.nn import functional as F
from torch.nn.init import xavier_normal
from torch.nn.modules.utils import _pair
from torch.utils.data import DataLoader
from tqdm import tqdm
import matplotlib.pyplot as plt
import scipy
import pandas as pd
import torchvision
from torchvision.transforms import Resize
import os
import pandas as pd
from PIL import Image
import torch
from torch.utils.data import Dataset
from torchvision.transforms import Resize
from numpy.fft import fft2, fftshift, fftfreq
We start by downloading one of the sessions available in the AllenSDK dataset.
data_dir = "./allendata"
#Again, we might want to use another folder if we do want to overwrite potential data we have
#data_dir = "./converted"
#From here on everything works with both AllenSDK and the MiniSDK!
manifest_path = os.path.join(data_dir, "manifest.json")
cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)
sessions = cache.get_session_table() #Returns a Pandas dataframe
sessions.head() #Sneak peek of the table
#Grab our (filtered) data from our favorite session. We take a female mice with nice unit count
session_id = 798911424
oursession = cache.get_session_data(session_id, timeout=3000)
In our analysis, we aim to investigate the neural response to visual stimuli, specifically focusing on natural scenes. Each natural scene is presented to the mouse 50 times, and during each presentation, the spiking times of a variable number of neurons are recorded simultaneously. Due to the inherent variability in neuronal activity from trial to trial, we have chosen to center our analysis on the average spike count of neurons in response to each image. This involves calculating the mean spike count across the multiple trials for each neuron.
We start by downloading the table of neurons' (unit_id) spiking times in response to one specific presentation of an image.
image_id = oursession.get_stimulus_table("natural_scenes").index.values
visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values
#Tables of the spikes happened during these presentations
visam_spikes = oursession.presentationwise_spike_times(
stimulus_presentation_ids=image_id,
unit_ids = visam_ids
)
visam_spikes
| stimulus_presentation_id | unit_id | time_since_stimulus_presentation_onset | |
|---|---|---|---|
| spike_time | |||
| 5909.794764 | 51355 | 951097547 | 0.000318 |
| 5909.796721 | 51355 | 951093080 | 0.002274 |
| 5909.801697 | 51355 | 951098928 | 0.007251 |
| 5909.803021 | 51355 | 951092303 | 0.008574 |
| 5909.805297 | 51355 | 951098560 | 0.010851 |
| ... | ... | ... | ... |
| 8569.511638 | 68228 | 951092303 | 0.244363 |
| 8569.512771 | 68228 | 951093688 | 0.245497 |
| 8569.515238 | 68228 | 951098021 | 0.247964 |
| 8569.516072 | 68228 | 951097682 | 0.248797 |
| 8569.517005 | 68228 | 951097869 | 0.249731 |
1184297 rows × 3 columns
We then download the table containing relevant information on the images of our dataset.
nat_scenes_temp = oursession.get_stimulus_table('natural_scenes')
#nat_scenes_temp.groupby('frame').count()
nat_scenes_temp
| stimulus_block | start_time | stop_time | frame | stimulus_name | duration | stimulus_condition_id | |
|---|---|---|---|---|---|---|---|
| stimulus_presentation_id | |||||||
| 51355 | 9.0 | 5909.794447 | 5910.044666 | 13.0 | natural_scenes | 0.250219 | 4908 |
| 51356 | 9.0 | 5910.044666 | 5910.294885 | 38.0 | natural_scenes | 0.250219 | 4909 |
| 51357 | 9.0 | 5910.294885 | 5910.545104 | 30.0 | natural_scenes | 0.250219 | 4910 |
| 51358 | 9.0 | 5910.545104 | 5910.795324 | 35.0 | natural_scenes | 0.250219 | 4911 |
| 51359 | 9.0 | 5910.795324 | 5911.045522 | 112.0 | natural_scenes | 0.250198 | 4912 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 68224 | 13.0 | 8568.266451 | 8568.516655 | 91.0 | natural_scenes | 0.250204 | 5005 |
| 68225 | 13.0 | 8568.516655 | 8568.766859 | 83.0 | natural_scenes | 0.250204 | 5022 |
| 68226 | 13.0 | 8568.766859 | 8569.017064 | 55.0 | natural_scenes | 0.250204 | 5014 |
| 68227 | 13.0 | 8569.017064 | 8569.267274 | 87.0 | natural_scenes | 0.250211 | 4955 |
| 68228 | 13.0 | 8569.267274 | 8569.517485 | 105.0 | natural_scenes | 0.250211 | 4942 |
5950 rows × 7 columns
We then perform some dataset manipulation to make it more readable, and we merge the two founded dataframes based on the id of the trial under consideration.
nat_scenes_temp = nat_scenes_temp.reset_index()
df = nat_scenes_temp.copy()
# Assuming df is your DataFrame with columns: 'stimuli', 'units', 'time'
# Sort the DataFrame by 'time' column
df_sorted = df.sort_values(by='start_time')
df_sorted.reset_index(inplace = True)
#final_data = df_first_presentation.merge(visam_spikes, on = 'stimulus_presentation_id', how = 'inner')
final_data = df_sorted.merge(visam_spikes, on = 'stimulus_presentation_id', how = 'inner')
final_data
| index | stimulus_presentation_id | stimulus_block | start_time | stop_time | x_position | phase | stimulus_name | size | y_position | contrast | temporal_frequency | spatial_frequency | orientation | duration | stimulus_condition_id | unit_id | time_since_stimulus_presentation_onset | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 1 | 0.0 | 84.942787 | 85.176306 | 20.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -10.0 | 0.8 | 4.0 | 0.08 | 90.0 | 0.233519 | 1 | 951093094 | 0.007928 |
| 1 | 0 | 1 | 0.0 | 84.942787 | 85.176306 | 20.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -10.0 | 0.8 | 4.0 | 0.08 | 90.0 | 0.233519 | 1 | 951098317 | 0.008801 |
| 2 | 0 | 1 | 0.0 | 84.942787 | 85.176306 | 20.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -10.0 | 0.8 | 4.0 | 0.08 | 90.0 | 0.233519 | 1 | 951092075 | 0.011728 |
| 3 | 0 | 1 | 0.0 | 84.942787 | 85.176306 | 20.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -10.0 | 0.8 | 4.0 | 0.08 | 90.0 | 0.233519 | 1 | 951092924 | 0.017895 |
| 4 | 0 | 1 | 0.0 | 84.942787 | 85.176306 | 20.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -10.0 | 0.8 | 4.0 | 0.08 | 90.0 | 0.233519 | 1 | 951098928 | 0.022368 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 425244 | 3644 | 3645 | 0.0 | 996.688034 | 996.938254 | -40.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -20.0 | 0.8 | 4.0 | 0.08 | 0.0 | 0.250221 | 123 | 951092973 | 0.225382 |
| 425245 | 3644 | 3645 | 0.0 | 996.688034 | 996.938254 | -40.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -20.0 | 0.8 | 4.0 | 0.08 | 0.0 | 0.250221 | 123 | 951092912 | 0.225849 |
| 425246 | 3644 | 3645 | 0.0 | 996.688034 | 996.938254 | -40.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -20.0 | 0.8 | 4.0 | 0.08 | 0.0 | 0.250221 | 123 | 951093795 | 0.237015 |
| 425247 | 3644 | 3645 | 0.0 | 996.688034 | 996.938254 | -40.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -20.0 | 0.8 | 4.0 | 0.08 | 0.0 | 0.250221 | 123 | 951093048 | 0.244382 |
| 425248 | 3644 | 3645 | 0.0 | 996.688034 | 996.938254 | -40.0 | [3644.93333333, 3644.93333333] | gabors | [20.0, 20.0] | -20.0 | 0.8 | 4.0 | 0.08 | 0.0 | 0.250221 | 123 | 951092303 | 0.249682 |
425249 rows × 18 columns
Then, we create the dataset with the spike count of neurons in response to every trial.
#Get the ids of the images and the units
#image_id = oursession.get_stimulus_table("natural_scenes").index.values
#visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values
image_id = final_data['stimulus_presentation_id'].values
visam_ids = final_data['unit_id'].values
print('done')
#Get the spikes during the session
#visam_spikes = oursession.presentationwise_spike_times(
# stimulus_presentation_ids=image_id,
# unit_ids = visam_ids
#)
visam_spikes = visam_spikes[visam_spikes['stimulus_presentation_id'].isin(image_id)]
visam_spikes = visam_spikes[visam_spikes['unit_id'].isin(visam_ids)]
print('done2')
#Add a new column to our table and fill it by counting number of rows with a stimulus presentation and unitid
visam_spikes["count"] = np.zeros(len(visam_spikes))
visam_spikes = visam_spikes.groupby(["stimulus_presentation_id", "unit_id"]).count()
#Employ "pivot table" to use the information from visam_spikes to generate a new table where
#we indicate index, columns, and fill it with the count.
spikes_presentation = pd.pivot_table(
visam_spikes,
values="time_since_stimulus_presentation_onset",
index="stimulus_presentation_id",
columns="unit_id",
fill_value=0.0,
aggfunc=np.sum
)
#This is 50 times faster!
spikes_presentation
done done2
| unit_id | 951092050 | 951092075 | 951092303 | 951092369 | 951092398 | 951092410 | 951092437 | 951092450 | 951092475 | 951092488 | ... | 951098773 | 951098807 | 951098850 | 951098871 | 951098928 | 951099272 | 951099306 | 951099320 | 951099475 | 951099491 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| stimulus_presentation_id | |||||||||||||||||||||
| 1 | 1 | 4 | 5 | 2 | 1 | 4 | 0 | 1 | 1 | 5 | ... | 8 | 1 | 1 | 3 | 4 | 0 | 2 | 0 | 0 | 0 |
| 2 | 1 | 2 | 1 | 2 | 0 | 4 | 0 | 0 | 2 | 1 | ... | 10 | 2 | 2 | 7 | 5 | 0 | 3 | 0 | 0 | 0 |
| 3 | 2 | 6 | 2 | 2 | 0 | 8 | 3 | 0 | 1 | 1 | ... | 5 | 1 | 0 | 6 | 5 | 0 | 2 | 1 | 0 | 0 |
| 4 | 1 | 3 | 2 | 3 | 1 | 1 | 2 | 1 | 0 | 7 | ... | 6 | 2 | 0 | 3 | 0 | 0 | 2 | 1 | 0 | 0 |
| 5 | 3 | 4 | 1 | 3 | 0 | 1 | 15 | 3 | 2 | 0 | ... | 6 | 4 | 1 | 5 | 3 | 0 | 0 | 0 | 0 | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3641 | 3 | 1 | 4 | 3 | 0 | 0 | 0 | 2 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 6 | 0 | 0 | 0 |
| 3642 | 3 | 1 | 3 | 2 | 4 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 2 | 0 | 5 | 0 | 0 | 0 |
| 3643 | 3 | 3 | 3 | 2 | 3 | 0 | 1 | 0 | 0 | 0 | ... | 1 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 |
| 3644 | 2 | 4 | 3 | 8 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 3 | 1 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 |
| 3645 | 2 | 2 | 5 | 2 | 1 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 2 | 0 | 0 |
3645 rows × 135 columns
We are now ready to average over a specific image. We firstly extract the image_id (frame column) from the original dataset, and merge it on the trial id. In this way we are then able to average over a specific frame id.
titles = final_data[['stimulus_presentation_id', 'frame']]
final_cnn = spikes_presentation.merge(titles, on = 'stimulus_presentation_id').drop_duplicates()
final_cnn.head()
| stimulus_presentation_id | 951092050 | 951092075 | 951092303 | 951092369 | 951092398 | 951092410 | 951092437 | 951092450 | 951092475 | ... | 951098807 | 951098850 | 951098871 | 951098928 | 951099272 | 951099306 | 951099320 | 951099475 | 951099491 | frame | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 51355 | 6 | 2 | 3 | 2 | 2 | 0 | 19 | 3 | 0 | ... | 9 | 1 | 0 | 11 | 0 | 0 | 4 | 0 | 1 | 13.0 |
| 372 | 51356 | 10 | 3 | 4 | 4 | 3 | 0 | 20 | 1 | 0 | ... | 12 | 0 | 2 | 8 | 0 | 0 | 2 | 0 | 5 | 38.0 |
| 728 | 51357 | 8 | 0 | 4 | 4 | 4 | 0 | 17 | 3 | 0 | ... | 9 | 1 | 0 | 8 | 0 | 0 | 0 | 0 | 3 | 30.0 |
| 1044 | 51358 | 9 | 5 | 2 | 5 | 8 | 1 | 10 | 9 | 0 | ... | 18 | 1 | 2 | 10 | 0 | 0 | 0 | 0 | 3 | 35.0 |
| 1444 | 51359 | 7 | 1 | 2 | 7 | 6 | 0 | 10 | 6 | 0 | ... | 10 | 0 | 1 | 10 | 0 | 0 | 4 | 0 | 6 | 112.0 |
5 rows × 137 columns
final_cnn = final_cnn.groupby('frame').mean()
final_cnn
| 951092050 | 951092075 | 951092303 | 951092369 | 951092398 | 951092410 | 951092437 | 951092450 | 951092475 | 951092488 | ... | 951098773 | 951098807 | 951098850 | 951098871 | 951098928 | 951099272 | 951099306 | 951099320 | 951099475 | 951099491 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| frame | |||||||||||||||||||||
| -1.0 | 3.48 | 2.06 | 4.54 | 5.28 | 0.52 | 0.68 | 2.84 | 1.22 | 0.16 | 3.42 | ... | 1.98 | 2.28 | 0.08 | 0.14 | 4.38 | 0.10 | 1.70 | 0.82 | 0.18 | 0.90 |
| 0.0 | 3.14 | 1.88 | 2.72 | 5.46 | 0.34 | 2.54 | 3.54 | 1.40 | 0.04 | 3.38 | ... | 1.62 | 2.60 | 0.12 | 0.04 | 3.80 | 0.16 | 0.80 | 0.46 | 0.18 | 0.74 |
| 1.0 | 3.14 | 2.18 | 3.88 | 5.04 | 0.72 | 1.02 | 4.20 | 1.04 | 0.04 | 2.26 | ... | 1.98 | 3.84 | 0.14 | 0.12 | 5.08 | 0.36 | 0.90 | 1.30 | 0.16 | 2.42 |
| 2.0 | 3.06 | 1.98 | 3.38 | 5.32 | 3.32 | 1.68 | 5.70 | 1.74 | 0.04 | 3.04 | ... | 2.24 | 3.14 | 0.10 | 0.18 | 6.20 | 0.14 | 1.12 | 0.78 | 0.08 | 1.08 |
| 3.0 | 3.16 | 1.94 | 3.46 | 4.34 | 1.20 | 2.26 | 5.28 | 1.44 | 0.08 | 2.48 | ... | 1.50 | 3.66 | 0.14 | 0.08 | 6.18 | 0.24 | 1.00 | 0.40 | 0.14 | 1.28 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 113.0 | 2.88 | 1.76 | 2.40 | 5.20 | 2.92 | 4.12 | 7.98 | 1.60 | 0.04 | 3.26 | ... | 4.20 | 5.30 | 0.30 | 0.36 | 6.86 | 0.10 | 3.98 | 1.38 | 0.02 | 0.46 |
| 114.0 | 3.14 | 2.66 | 3.08 | 4.64 | 0.90 | 0.62 | 4.62 | 1.14 | 0.06 | 1.98 | ... | 2.18 | 3.76 | 0.08 | 0.18 | 5.44 | 0.22 | 1.32 | 0.82 | 0.20 | 0.80 |
| 115.0 | 3.68 | 2.02 | 3.30 | 5.28 | 0.68 | 2.66 | 6.52 | 1.36 | 0.12 | 2.54 | ... | 2.12 | 3.70 | 0.20 | 0.14 | 6.00 | 0.16 | 2.74 | 0.64 | 0.14 | 1.00 |
| 116.0 | 4.00 | 2.92 | 4.00 | 5.52 | 0.58 | 1.30 | 4.98 | 1.40 | 0.12 | 2.30 | ... | 2.52 | 2.94 | 0.10 | 0.16 | 4.44 | 0.14 | 0.78 | 0.44 | 0.06 | 1.16 |
| 117.0 | 3.54 | 2.50 | 3.14 | 5.16 | 0.90 | 1.94 | 7.12 | 1.32 | 0.06 | 2.04 | ... | 2.68 | 4.16 | 0.02 | 0.10 | 5.48 | 0.12 | 1.28 | 1.86 | 0.24 | 0.98 |
119 rows × 135 columns
We save the dataset for later uses.
final_cnn.to_csv('dataset_for_agne_mean.csv')
To train the Convolutional Neural Network (CNN) and identify the Most Excitatory Inputs (MEIs), we constructed a spiking activity dataset.
This dataset is organized in a matrix format, where columns represent unique neuron IDs, rows represent unique image IDs, and each cell contains the average firing rate (spikes per unit time) of a specific neuron in response to a specific image across all trials.
df = pd.read_csv("dataset_for_agne_mean.csv")
df
| frame | stimulus_presentation_id | 951092050 | 951092075 | 951092303 | 951092369 | 951092398 | 951092410 | 951092437 | 951092450 | ... | 951098773 | 951098807 | 951098850 | 951098871 | 951098928 | 951099272 | 951099306 | 951099320 | 951099475 | 951099491 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -1.0 | 60469.34 | 3.48 | 2.06 | 4.54 | 5.28 | 0.52 | 0.68 | 2.84 | 1.22 | ... | 1.98 | 2.28 | 0.08 | 0.14 | 4.38 | 0.10 | 1.70 | 0.82 | 0.18 | 0.90 |
| 1 | 0.0 | 58409.40 | 3.14 | 1.88 | 2.72 | 5.46 | 0.34 | 2.54 | 3.54 | 1.40 | ... | 1.62 | 2.60 | 0.12 | 0.04 | 3.80 | 0.16 | 0.80 | 0.46 | 0.18 | 0.74 |
| 2 | 1.0 | 59487.90 | 3.14 | 2.18 | 3.88 | 5.04 | 0.72 | 1.02 | 4.20 | 1.04 | ... | 1.98 | 3.84 | 0.14 | 0.12 | 5.08 | 0.36 | 0.90 | 1.30 | 0.16 | 2.42 |
| 3 | 2.0 | 57120.46 | 3.06 | 1.98 | 3.38 | 5.32 | 3.32 | 1.68 | 5.70 | 1.74 | ... | 2.24 | 3.14 | 0.10 | 0.18 | 6.20 | 0.14 | 1.12 | 0.78 | 0.08 | 1.08 |
| 4 | 3.0 | 59864.54 | 3.16 | 1.94 | 3.46 | 4.34 | 1.20 | 2.26 | 5.28 | 1.44 | ... | 1.50 | 3.66 | 0.14 | 0.08 | 6.18 | 0.24 | 1.00 | 0.40 | 0.14 | 1.28 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 114 | 113.0 | 59083.06 | 2.88 | 1.76 | 2.40 | 5.20 | 2.92 | 4.12 | 7.98 | 1.60 | ... | 4.20 | 5.30 | 0.30 | 0.36 | 6.86 | 0.10 | 3.98 | 1.38 | 0.02 | 0.46 |
| 115 | 114.0 | 57132.22 | 3.14 | 2.66 | 3.08 | 4.64 | 0.90 | 0.62 | 4.62 | 1.14 | ... | 2.18 | 3.76 | 0.08 | 0.18 | 5.44 | 0.22 | 1.32 | 0.82 | 0.20 | 0.80 |
| 116 | 115.0 | 59239.70 | 3.68 | 2.02 | 3.30 | 5.28 | 0.68 | 2.66 | 6.52 | 1.36 | ... | 2.12 | 3.70 | 0.20 | 0.14 | 6.00 | 0.16 | 2.74 | 0.64 | 0.14 | 1.00 |
| 117 | 116.0 | 58072.78 | 4.00 | 2.92 | 4.00 | 5.52 | 0.58 | 1.30 | 4.98 | 1.40 | ... | 2.52 | 2.94 | 0.10 | 0.16 | 4.44 | 0.14 | 0.78 | 0.44 | 0.06 | 1.16 |
| 118 | 117.0 | 57424.92 | 3.54 | 2.50 | 3.14 | 5.16 | 0.90 | 1.94 | 7.12 | 1.32 | ... | 2.68 | 4.16 | 0.02 | 0.10 | 5.48 | 0.12 | 1.28 | 1.86 | 0.24 | 0.98 |
119 rows × 137 columns
data_dir = "./allendata"
manifest_path = os.path.join(data_dir, "manifest.json")
cache = EcephysProjectCache.from_warehouse(manifest=manifest_path)
sessions = cache.get_session_table() #Returns a Pandas dataframe
session_id = 798911424
oursession = cache.get_session_data(session_id, timeout=3000)
nat_scenes_temp = oursession.get_stimulus_table('natural_scenes')
nat_scenes_temp = nat_scenes_temp.reset_index()
nat_scenes_temp
| stimulus_presentation_id | stimulus_block | start_time | stop_time | stimulus_name | frame | duration | stimulus_condition_id | |
|---|---|---|---|---|---|---|---|---|
| 0 | 51355 | 9.0 | 5909.794447 | 5910.044666 | natural_scenes | 13.0 | 0.250219 | 4908 |
| 1 | 51356 | 9.0 | 5910.044666 | 5910.294885 | natural_scenes | 38.0 | 0.250219 | 4909 |
| 2 | 51357 | 9.0 | 5910.294885 | 5910.545104 | natural_scenes | 30.0 | 0.250219 | 4910 |
| 3 | 51358 | 9.0 | 5910.545104 | 5910.795324 | natural_scenes | 35.0 | 0.250219 | 4911 |
| 4 | 51359 | 9.0 | 5910.795324 | 5911.045522 | natural_scenes | 112.0 | 0.250198 | 4912 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 5945 | 68224 | 13.0 | 8568.266451 | 8568.516655 | natural_scenes | 91.0 | 0.250204 | 5005 |
| 5946 | 68225 | 13.0 | 8568.516655 | 8568.766859 | natural_scenes | 83.0 | 0.250204 | 5022 |
| 5947 | 68226 | 13.0 | 8568.766859 | 8569.017064 | natural_scenes | 55.0 | 0.250204 | 5014 |
| 5948 | 68227 | 13.0 | 8569.017064 | 8569.267274 | natural_scenes | 87.0 | 0.250211 | 4955 |
| 5949 | 68228 | 13.0 | 8569.267274 | 8569.517485 | natural_scenes | 105.0 | 0.250211 | 4942 |
5950 rows × 8 columns
image_id = oursession.get_stimulus_table("natural_scenes").index.values
visam_ids = oursession.units[oursession.units["ecephys_structure_acronym"]=="VISam"].index.values
visam_spikes = oursession.presentationwise_spike_times(
stimulus_presentation_ids=image_id,
unit_ids = visam_ids
)
visam_spikes.reset_index(inplace=True)
visam_spikes
| spike_time | stimulus_presentation_id | unit_id | time_since_stimulus_presentation_onset | |
|---|---|---|---|---|
| 0 | 5909.794764 | 51355 | 951097547 | 0.000318 |
| 1 | 5909.796721 | 51355 | 951093080 | 0.002274 |
| 2 | 5909.801697 | 51355 | 951098928 | 0.007251 |
| 3 | 5909.803021 | 51355 | 951092303 | 0.008574 |
| 4 | 5909.805297 | 51355 | 951098560 | 0.010851 |
| ... | ... | ... | ... | ... |
| 1184292 | 8569.511638 | 68228 | 951092303 | 0.244363 |
| 1184293 | 8569.512771 | 68228 | 951093688 | 0.245497 |
| 1184294 | 8569.515238 | 68228 | 951098021 | 0.247964 |
| 1184295 | 8569.516072 | 68228 | 951097682 | 0.248797 |
| 1184296 | 8569.517005 | 68228 | 951097869 | 0.249731 |
1184297 rows × 4 columns
To identify the Most Important Neurons (MINs) using the oracle correlation method, we constructed a separate dataset.
This dataset contains trial IDs, unique neuron IDs, unique image IDs, and the corresponding spike count for each specific neuron in response to a particular image during a specific trial.
# We're interested in the spikes that are in response to the natural scenes
df_new = visam_spikes.copy()
# Drop the columns that we don't need
df_new.drop(columns=["time_since_stimulus_presentation_onset"], inplace=True)
# Count the number of spikes per stimulus presentation
df_new = df_new.groupby(['stimulus_presentation_id', 'unit_id']).count()
df_new.reset_index(inplace=True)
df_new
| stimulus_presentation_id | unit_id | spike_time | |
|---|---|---|---|
| 0 | 51355 | 951092050 | 6 |
| 1 | 51355 | 951092075 | 2 |
| 2 | 51355 | 951092303 | 3 |
| 3 | 51355 | 951092369 | 2 |
| 4 | 51355 | 951092398 | 2 |
| ... | ... | ... | ... |
| 386344 | 68228 | 951098773 | 3 |
| 386345 | 68228 | 951098807 | 2 |
| 386346 | 68228 | 951098850 | 1 |
| 386347 | 68228 | 951098928 | 7 |
| 386348 | 68228 | 951099320 | 6 |
386349 rows × 3 columns
# Take natural scenes in order to take also the image IDs
df_natural = oursession.get_stimulus_table('natural_scenes')
df_natural.reset_index(inplace = True)
# Drop useless columns
df_natural.drop(columns=['stimulus_block', 'start_time', 'stop_time', 'stimulus_name', 'duration', 'stimulus_condition_id'], inplace=True)
# Rename the columns
df_natural.rename(columns={'stimulus_presentation_id': 'stimulus_presentation_id', 'frame':'image_id'}, inplace=True)
df_natural
| stimulus_presentation_id | image_id | |
|---|---|---|
| 0 | 51355 | 13.0 |
| 1 | 51356 | 38.0 |
| 2 | 51357 | 30.0 |
| 3 | 51358 | 35.0 |
| 4 | 51359 | 112.0 |
| ... | ... | ... |
| 5945 | 68224 | 91.0 |
| 5946 | 68225 | 83.0 |
| 5947 | 68226 | 55.0 |
| 5948 | 68227 | 87.0 |
| 5949 | 68228 | 105.0 |
5950 rows × 2 columns
# Merge the two datasets
final_df = df_new.merge(df_natural, how ='left', on='stimulus_presentation_id' )
final_df
| stimulus_presentation_id | unit_id | spike_time | image_id | |
|---|---|---|---|---|
| 0 | 51355 | 951092050 | 6 | 13.0 |
| 1 | 51355 | 951092075 | 2 | 13.0 |
| 2 | 51355 | 951092303 | 3 | 13.0 |
| 3 | 51355 | 951092369 | 2 | 13.0 |
| 4 | 51355 | 951092398 | 2 | 13.0 |
| ... | ... | ... | ... | ... |
| 386344 | 68228 | 951098773 | 3 | 105.0 |
| 386345 | 68228 | 951098807 | 2 | 105.0 |
| 386346 | 68228 | 951098850 | 1 | 105.0 |
| 386347 | 68228 | 951098928 | 7 | 105.0 |
| 386348 | 68228 | 951099320 | 6 | 105.0 |
386349 rows × 4 columns
# Rename the columns
final_df.rename(columns={'stimulus_presentation_id': "trial_id", "unit_id": "neuron_id", "spike_time": "spike_count", "image_id": "image_id"}, inplace = True)
# Convert image_id to integer
convert_dict = {'image_id':int}
final_df = final_df.astype(convert_dict)
final_df
| trial_id | neuron_id | spike_count | image_id | |
|---|---|---|---|---|
| 0 | 51355 | 951092050 | 6 | 13 |
| 1 | 51355 | 951092075 | 2 | 13 |
| 2 | 51355 | 951092303 | 3 | 13 |
| 3 | 51355 | 951092369 | 2 | 13 |
| 4 | 51355 | 951092398 | 2 | 13 |
| ... | ... | ... | ... | ... |
| 386344 | 68228 | 951098773 | 3 | 105 |
| 386345 | 68228 | 951098807 | 2 | 105 |
| 386346 | 68228 | 951098850 | 1 | 105 |
| 386347 | 68228 | 951098928 | 7 | 105 |
| 386348 | 68228 | 951099320 | 6 | 105 |
386349 rows × 4 columns
final_df.to_csv("oracle_df.csv", index=False)
min_time = visam_spikes['spike_time'].min()
max_time = visam_spikes['spike_time'].max()
binsize = 0.1 # BINSIZE (!!!!)
bins = np.arange(min_time, max_time, binsize)
min_time, max_time
(5909.794764092189, 8569.517004859086)
def count_spikes(group):
counts, _ = np.histogram(group['spike_time'], bins=bins)
return counts
spike_counts = visam_spikes.groupby('unit_id').apply(count_spikes)
# Convert the Series of arrays into a 2D NumPy array
spike_counts_matrix = np.stack(spike_counts.values)
spike_counts_matrix
array([[2, 3, 2, ..., 2, 1, 2],
[1, 1, 2, ..., 0, 0, 1],
[1, 1, 3, ..., 2, 2, 2],
...,
[1, 1, 2, ..., 1, 3, 3],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 2, ..., 0, 0, 0]])
spike_counts_zscored = zscore(spike_counts_matrix, axis=1)
spike_counts_zscored.shape
(135, 26597)
model = Rastermap(n_clusters=None, # None turns off clustering and sorts single neurons
n_PCs=64, # use fewer PCs than neurons
locality=0.15, # some locality in sorting (this is a value from 0-1)
time_lag_window=15, # use future timepoints to compute correlation
grid_upsample=0, # 0 turns off upsampling since we're using single neurons
).fit(spike_counts_zscored)
y = model.embedding # neurons x 1
isort = model.isort
normalizing data across axis=1 projecting out mean along axis=0 data normalized, 0.06sec sorting activity: 135 valid samples by 26597 timepoints n_PCs = 64 computed, 4.75sec skipping clustering, n_clusters is None clusters sorted, time 135.09sec rastermap complete, time 135.10sec
isort = model.isort
embedding = model.embedding
xmin = 22000
xmax = 28000
fig, ax = plt.subplots(figsize=(12, 8), dpi=200)
im = ax.imshow(spike_counts_zscored[isort, xmin:xmax], cmap="gray_r", vmin=0, vmax=1.2, aspect="auto")
ax.set_xlabel("Time (bins)")
ax.set_ylabel("Neurons (sorted)")
ax.set_title("Raster Plot of Neural Activity")
cbar = plt.colorbar(im, ax=ax, pad=0.01)
cbar.set_label('Spike count')
plt.show()
The mouse has been shown natural images three different times, in this plot we focused on the last display.
We can now focus on the single blocks
# Define different xmin and xmax for each plot
xmin = [0, 8000, 22000]
xmax = [4800, 12000, 28000]
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 6), dpi=200)
for i, ax in enumerate(axes.flat):
# Plot each segment of the data
im = ax.imshow(spike_counts_zscored[isort, xmin[i]:xmax[i]], cmap="gray_r", vmin=0, vmax=1.2, aspect="auto")
ax.set_xlabel("Time (bins)")
ax.set_ylabel("Neurons (sorted)")
ax.set_title(f"Raster Plot from {xmin[i]} to {xmax[i]} bins")
plt.tight_layout()
plt.show()
This plot shows that there are indeed some neurons which present a pattern behaviour in neural response, which indicates that there are good possibilities to select a neuron which will work well with our CNN.
However let us note that there are neurons whose behaviour seems more noisy, hence selecting a good neurons is of utmost importance.
Importing the dataset, dropping useless data, transposing in order to have neurons as rows, columns as images
data = pd.read_csv('dataset_for_agne_mean.csv')
data.drop(['frame', 'stimulus_presentation_id'], axis=1, inplace=True)
data = data.transpose()
Plotting the correlation matrix as a heatmap (transposing again in order to do it on neurons)
correlation_matrix = data.transpose().corr()
# Plotting the correlation matrix as a heatmap
plt.figure(figsize=(10, 10)) # Adjust the size as needed
sns.heatmap(correlation_matrix, cmap='coolwarm', center=0)
plt.show()
No particular ambiguous correlation can be seen from this plot.
Now let's try to find for each neuron, the most correlated and the least correlated neuron.
In order to avoid problems, we fill the diagonal with NaN values, since for each neuron, the correlation with itself is always $1.0$
np.fill_diagonal(correlation_matrix.values, np.nan)
highest_correlations = correlation_matrix.max()
lowest_correlation = correlation_matrix.min()
most_correlated_neurons = correlation_matrix.idxmax()
least_correlated_neurons = correlation_matrix.idxmin()
results = pd.DataFrame({
'Neuron': highest_correlations.index,
'Highest Correlation': highest_correlations.values,
'Most Correlated Neuron': most_correlated_neurons.values,
'Lowest Correlation': lowest_correlation.values,
'Least Correlated Neuron': least_correlated_neurons.values
})
print(results)
Neuron Highest Correlation Most Correlated Neuron \
0 951092050 0.297214 951093080
1 951092075 0.309912 951098279
2 951092303 0.416028 951093160
3 951092369 0.256688 951098530
4 951092398 0.655551 951092410
.. ... ... ...
130 951099272 0.342826 951098414
131 951099306 0.797598 951094081
132 951099320 0.462756 951098242
133 951099475 0.392041 951098430
134 951099491 0.403545 951092303
Lowest Correlation Least Correlated Neuron
0 -0.242279 951092488
1 -0.239617 951093510
2 -0.476467 951092949
3 -0.240787 951093608
4 -0.418246 951092303
.. ... ...
130 -0.307064 951092890
131 -0.576034 951097947
132 -0.462980 951097826
133 -0.257062 951098075
134 -0.369061 951097452
[135 rows x 5 columns]
Let's see the highest correlated neurons by fixing a threshold of $0.7$.
highly_correlated_neurons = results[results['Highest Correlation'] > 0.7]
highly_correlated_neurons
| Neuron | Highest Correlation | Most Correlated Neuron | Lowest Correlation | Least Correlated Neuron | |
|---|---|---|---|---|---|
| 32 | 951093125 | 0.905601 | 951098242 | -0.553941 | 951097947 |
| 47 | 951093608 | 0.736594 | 951093821 | -0.401066 | 951098279 |
| 56 | 951093821 | 0.736594 | 951093608 | -0.429395 | 951092303 |
| 60 | 951093928 | 0.750732 | 951098242 | -0.526685 | 951097947 |
| 64 | 951094081 | 0.848041 | 951098242 | -0.548484 | 951097947 |
| 94 | 951097826 | 0.703051 | 951098560 | -0.462980 | 951099320 |
| 105 | 951098242 | 0.905601 | 951093125 | -0.654761 | 951097947 |
| 123 | 951098560 | 0.703051 | 951097826 | -0.372543 | 951098303 |
| 131 | 951099306 | 0.797598 | 951094081 | -0.576034 | 951097947 |
Let's see the lowest correlated neurons (i.e., negatively correlated)
lowly_correlated_neurons = results[results['Lowest Correlation'] < -0.55]
lowly_correlated_neurons
| Neuron | Highest Correlation | Most Correlated Neuron | Lowest Correlation | Least Correlated Neuron | |
|---|---|---|---|---|---|
| 32 | 951093125 | 0.905601 | 951098242 | -0.553941 | 951097947 |
| 99 | 951097947 | 0.368022 | 951093608 | -0.654761 | 951098242 |
| 105 | 951098242 | 0.905601 | 951093125 | -0.654761 | 951097947 |
| 131 | 951099306 | 0.797598 | 951094081 | -0.576034 | 951097947 |
This result is weird, since by fixing a threshold of $-0.55$, we see that $4$ neurons are more negatively correlated than that, with the highest being $-0.654761$.
df = data.copy()
iso_forest = IsolationForest(n_estimators=100, contamination='auto', random_state=42)
iso_forest.fit(df.values)
anomalies = iso_forest.predict(df.values)
df['anomaly_score'] = iso_forest.decision_function(df.values)
df['is_anomaly'] = anomalies
anomalous_neurons = df[df['is_anomaly'] == -1]
neuron_entropies = data.apply(lambda x: entropy(x.value_counts()), axis=1)
# Dataset with anomaly values, entropy, correlation, and neuron with the most correlation
df_entropy_anomaly = pd.DataFrame({
'Entropy' : neuron_entropies,
'Anomaly' : df['is_anomaly'],
'Highest reported correlation with another neuron': highest_correlations.values,
'Most Correlated Neuron': most_correlated_neurons.values,
'Variance': df.var(axis=1)
})
df_entropy_anomaly.sort_values(by='Entropy', ascending=False, inplace=True)
df_entropy_anomaly.reset_index(inplace=True)
type_dict = {'index': int}
df_entropy_anomaly = df_entropy_anomaly.astype(type_dict)
df_entropy_anomaly
| index | Entropy | Anomaly | Highest reported correlation with another neuron | Most Correlated Neuron | Variance | |
|---|---|---|---|---|---|---|
| 0 | 951093821 | 4.465900 | -1 | 0.736594 | 951093608 | 2.694286 |
| 1 | 951093608 | 4.435348 | -1 | 0.736594 | 951093821 | 2.418180 |
| 2 | 951092437 | 4.423117 | -1 | 0.515452 | 951097547 | 1.912515 |
| 3 | 951093283 | 4.400400 | -1 | 0.414326 | 951092437 | 1.576702 |
| 4 | 951098928 | 4.362596 | -1 | 0.503125 | 951092437 | 1.434995 |
| ... | ... | ... | ... | ... | ... | ... |
| 130 | 951093810 | 2.104046 | 1 | 0.384974 | 951094050 | 0.011959 |
| 131 | 951093510 | 2.061717 | 1 | 0.344498 | 951094573 | 0.012383 |
| 132 | 951098317 | 2.045806 | 1 | 0.299732 | 951093510 | 0.009211 |
| 133 | 951098372 | 2.014949 | 1 | 0.397437 | 951093283 | 0.009025 |
| 134 | 951094036 | 1.533163 | 1 | 0.339248 | 951099306 | 0.014631 |
135 rows × 6 columns
df_entropy_anomaly['Anomaly'].value_counts()
1 116 -1 19 Name: Anomaly, dtype: int64
df_entropy_anomaly.head(20)[df_entropy_anomaly['Anomaly']==-1].shape
/var/folders/_2/p0xm13rx0yq8v664xb_zz8tc0000gn/T/ipykernel_78987/1485505297.py:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index. df_entropy_anomaly.head(20)[df_entropy_anomaly['Anomaly']==-1].shape
(14, 6)
df_entropy_anomaly.sort_values(by='Variance', ascending=False).head(20)[df_entropy_anomaly['Anomaly']==-1].shape
/var/folders/_2/p0xm13rx0yq8v664xb_zz8tc0000gn/T/ipykernel_78987/2583719294.py:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index. df_entropy_anomaly.sort_values(by='Variance', ascending=False).head(20)[df_entropy_anomaly['Anomaly']==-1].shape
(15, 6)
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly, x='Entropy', y='Variance', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Entropy vs Variance of Neuron Responses')
plt.xlabel('Entropy')
plt.ylabel('Variance')
plt.legend(title='Anomaly')
plt.show()
It's interesting to note that the isolation forest identified as anomalies neurons with high variance and high entropy.
df_oracle = pd.read_csv('oracle_df.csv')
df_oracle
| trial_id | neuron_id | spike_count | image_id | |
|---|---|---|---|---|
| 0 | 51355 | 951092050 | 6 | 13 |
| 1 | 51355 | 951092075 | 2 | 13 |
| 2 | 51355 | 951092303 | 3 | 13 |
| 3 | 51355 | 951092369 | 2 | 13 |
| 4 | 51355 | 951092398 | 2 | 13 |
| ... | ... | ... | ... | ... |
| 386344 | 68228 | 951098773 | 3 | 105 |
| 386345 | 68228 | 951098807 | 2 | 105 |
| 386346 | 68228 | 951098850 | 1 | 105 |
| 386347 | 68228 | 951098928 | 7 | 105 |
| 386348 | 68228 | 951099320 | 6 | 105 |
386349 rows × 4 columns
def super_oracle(df):
corr = dict()
# Total number of images
total_images = len(df['image_id'].unique())
# Take a neuron
for neuron_id in df['neuron_id'].unique():
responses_list = []
loom_list = []
# Image counter
image_counter = 0
# Take an image which is shown to the neuron
for image_id in df[(df['neuron_id'] == neuron_id)]['image_id'].unique():
# Skip if there is only one trial for the image
if len(df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id)]['trial_id'].unique()) == 1:
continue
# Take a random trial for the image and neuron combination
random_trial = random.choice(df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id)]['trial_id'].unique())
# Take the neuron response to the random trial
response = df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id) & (df['trial_id']==random_trial)]['spike_count'].values[0]
# Compute the mean over all the other trials for the image
loom_image = df[(df['neuron_id'] == neuron_id) & (df['image_id'] == image_id) & (df['trial_id']!=random_trial)]['spike_count'].values.mean()
# Increment the image counter
image_counter += 1
# Append
responses_list.append(response)
loom_list.append(loom_image)
# Compute the correlation between the leave-one-out-mean and the actual response for the neuron on that image
try:
correlation = pearsonr(responses_list, loom_list)[0]
except:
print(f"Neuron {neuron_id} has a problem")
return neuron_id
# Compute the image fraction
image_fraction = image_counter / total_images
corr[neuron_id] = correlation*image_fraction
return pd.DataFrame.from_dict(corr, orient='index', columns=['super_oracle_correlation'])
cor_super = super_oracle(df_oracle)
cor_super
| super_oracle_correlation | |
|---|---|
| 951092050 | -0.115015 |
| 951092075 | -0.016209 |
| 951092303 | 0.134748 |
| 951092369 | 0.089028 |
| 951092398 | 0.517139 |
| ... | ... |
| 951099272 | 0.037873 |
| 951099475 | -0.004608 |
| 951098292 | 0.133259 |
| 951094453 | 0.268118 |
| 951098317 | -0.029089 |
135 rows × 1 columns
cor_super.reset_index(inplace=True)
cor_super.sort_values(by='super_oracle_correlation', ascending=False, inplace=True)
cor_super
| index | super_oracle_correlation | |
|---|---|---|
| 41 | 951094050 | 0.651761 |
| 108 | 951093125 | 0.579137 |
| 102 | 951098242 | 0.578967 |
| 30 | 951093608 | 0.569689 |
| 4 | 951092398 | 0.517139 |
| ... | ... | ... |
| 7 | 951092488 | -0.093438 |
| 84 | 951093416 | -0.096325 |
| 106 | 951093020 | -0.105021 |
| 0 | 951092050 | -0.115015 |
| 28 | 951093565 | -0.163587 |
135 rows × 2 columns
def average_super_oracle(df, iterations=10):
rank_sums = {} # Stores sum of ranks
correlation_sums = {} # Stores sum of correlations
total_neuron_ids = df['neuron_id'].unique()
# Initialize dictionaries
for neuron_id in total_neuron_ids:
rank_sums[neuron_id] = 0
correlation_sums[neuron_id] = 0
for _ in range(iterations):
# Run the super_oracle function
result_df = super_oracle(df)
if type(result_df) != pd.DataFrame:
print(f"Error in neuron {result_df}")
return 0
# Sort results by correlation in descending order and reset index
sorted_df = result_df.sort_values(by='super_oracle_correlation', ascending=False).reset_index()
# Record ranks and sum up correlations
for index, row in sorted_df.iterrows():
neuron_id = row['index']
correlation = row['super_oracle_correlation']
rank_sums[neuron_id] += index # index is the rank
correlation_sums[neuron_id] += correlation
# Prepare to return the results
average_ranks = {neuron_id: rank_sums[neuron_id] / iterations for neuron_id in total_neuron_ids}
average_correlations = {neuron_id: correlation_sums[neuron_id] / iterations for neuron_id in total_neuron_ids}
# Create DataFrame from the computed averages
summary_df = pd.DataFrame.from_dict({
'Average Rank': average_ranks,
'Average Correlation': average_correlations
})
# Sort by average rank
summary_df = summary_df.sort_values(by='Average Rank')
return summary_df
av = average_super_oracle(df_oracle, iterations=50)
av
| Average Rank | Average Correlation | |
|---|---|---|
| 951098242 | 1.44 | 0.602705 |
| 951093125 | 2.94 | 0.571175 |
| 951092798 | 5.70 | 0.489572 |
| 951093846 | 6.88 | 0.499745 |
| 951093821 | 7.84 | 0.452857 |
| ... | ... | ... |
| 951098317 | 116.82 | -0.043106 |
| 951093858 | 116.84 | -0.044025 |
| 951093565 | 117.58 | -0.081306 |
| 951093416 | 122.58 | -0.082603 |
| 951092673 | 126.20 | -0.085904 |
135 rows × 2 columns
cor_super = av.copy()
cor_super.drop(columns=['Average Rank'], inplace=True)
cor_super.rename(columns={'Average Correlation': 'oracle_correlation'}, inplace=True)
cor_super.reset_index(inplace=True)
cor_super
| index | oracle_correlation | |
|---|---|---|
| 0 | 951098242 | 0.602705 |
| 1 | 951093125 | 0.571175 |
| 2 | 951092798 | 0.489572 |
| 3 | 951093846 | 0.499745 |
| 4 | 951093821 | 0.452857 |
| ... | ... | ... |
| 130 | 951098317 | -0.043106 |
| 131 | 951093858 | -0.044025 |
| 132 | 951093565 | -0.081306 |
| 133 | 951093416 | -0.082603 |
| 134 | 951092673 | -0.085904 |
135 rows × 2 columns
df_entropy_anomaly_oracle = pd.merge(df_entropy_anomaly, cor_super, left_on='index', right_on='index')
df_entropy_anomaly_oracle.sort_values(by='oracle_correlation', ascending=False, inplace=True)
df_entropy_anomaly_oracle.reset_index(inplace=True)
df_entropy_anomaly_oracle.drop(['level_0'], axis=1, inplace=True)
df_entropy_anomaly_oracle
| index | Entropy | Anomaly | Highest reported correlation with another neuron | Most Correlated Neuron | Variance | oracle_correlation | |
|---|---|---|---|---|---|---|---|
| 0 | 951098242 | 3.992573 | 1 | 0.905601 | 951093125 | 0.968955 | 0.602705 |
| 1 | 951093125 | 3.716815 | 1 | 0.905601 | 951098242 | 0.497883 | 0.571175 |
| 2 | 951093846 | 3.899434 | 1 | 0.677264 | 951092450 | 0.769341 | 0.499745 |
| 3 | 951092798 | 4.129662 | 1 | 0.482546 | 951092949 | 0.725587 | 0.489572 |
| 4 | 951093821 | 4.465900 | -1 | 0.736594 | 951093608 | 2.694286 | 0.452857 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 130 | 951093858 | 2.302911 | 1 | 0.368254 | 951094654 | 0.011904 | -0.044025 |
| 131 | 951093565 | 3.004287 | 1 | 0.279689 | 951094050 | 0.016774 | -0.081306 |
| 132 | 951093416 | 3.134690 | 1 | 0.436531 | 951099306 | 0.027655 | -0.082603 |
| 133 | 951092673 | 2.696601 | 1 | 0.410253 | 951097530 | 0.019744 | -0.085904 |
| 134 | 951098372 | 2.014949 | 1 | 0.397437 | 951093283 | 0.009025 | NaN |
135 rows × 7 columns
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly_oracle[df_entropy_anomaly_oracle['oracle_correlation']>0], x='Entropy', y='oracle_correlation', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Entropy vs Oracle Correlation of Neuron Responses')
plt.xlabel('Entropy')
plt.ylabel('Oracle Correlation')
plt.legend(title='Anomaly')
plt.show()
plt.figure(figsize=(10, 6))
sns.scatterplot(data=df_entropy_anomaly_oracle, x='Variance', y='oracle_correlation', hue='Anomaly', style='Anomaly', alpha=0.7)
plt.title('Variance vs Oracle Correlation of Neuron Responses')
plt.xlabel('Variance')
plt.ylabel('Oracle Correlation')
plt.legend(title='Anomaly')
plt.show()
df_entropy_anomaly_oracle.sort_values(by='oracle_correlation', ascending= False)
| index | Entropy | Anomaly | Highest reported correlation with another neuron | Most Correlated Neuron | Variance | oracle_correlation | Color | |
|---|---|---|---|---|---|---|---|---|
| 0 | 951098242 | 3.992573 | 1 | 0.905601 | 951093125 | 0.968955 | 0.602705 | Other Neurons |
| 1 | 951093125 | 3.716815 | 1 | 0.905601 | 951098242 | 0.497883 | 0.571175 | Other Neurons |
| 2 | 951093846 | 3.899434 | 1 | 0.677264 | 951092450 | 0.769341 | 0.499745 | Other Neurons |
| 3 | 951092798 | 4.129662 | 1 | 0.482546 | 951092949 | 0.725587 | 0.489572 | Other Neurons |
| 4 | 951093821 | 4.465900 | -1 | 0.736594 | 951093608 | 2.694286 | 0.452857 | Other Neurons |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 130 | 951093858 | 2.302911 | 1 | 0.368254 | 951094654 | 0.011904 | -0.044025 | Other Neurons |
| 131 | 951093565 | 3.004287 | 1 | 0.279689 | 951094050 | 0.016774 | -0.081306 | Other Neurons |
| 132 | 951093416 | 3.134690 | 1 | 0.436531 | 951099306 | 0.027655 | -0.082603 | Other Neurons |
| 133 | 951092673 | 2.696601 | 1 | 0.410253 | 951097530 | 0.019744 | -0.085904 | Other Neurons |
| 134 | 951098372 | 2.014949 | 1 | 0.397437 | 951093283 | 0.009025 | NaN | Other Neurons |
135 rows × 8 columns
# Take neurons in the 60th percentile of the oracle correlation and 50th percentile of the variance
top_oracle_neurons = df_entropy_anomaly_oracle[(df_entropy_anomaly_oracle['oracle_correlation'] > df_entropy_anomaly_oracle['oracle_correlation'].quantile(0.6)) & (df_entropy_anomaly_oracle['Variance'] > df_entropy_anomaly_oracle['Variance'].quantile(0.5))]
top_oracle_neurons
| index | Entropy | Anomaly | Highest reported correlation with another neuron | Most Correlated Neuron | Variance | oracle_correlation | Color | |
|---|---|---|---|---|---|---|---|---|
| 0 | 951098242 | 3.992573 | 1 | 0.905601 | 951093125 | 0.968955 | 0.602705 | Other Neurons |
| 1 | 951093125 | 3.716815 | 1 | 0.905601 | 951098242 | 0.497883 | 0.571175 | Other Neurons |
| 2 | 951093846 | 3.899434 | 1 | 0.677264 | 951092450 | 0.769341 | 0.499745 | Other Neurons |
| 3 | 951092798 | 4.129662 | 1 | 0.482546 | 951092949 | 0.725587 | 0.489572 | Other Neurons |
| 4 | 951093821 | 4.465900 | -1 | 0.736594 | 951093608 | 2.694286 | 0.452857 | Other Neurons |
| 5 | 951097947 | 4.232909 | -1 | 0.368022 | 951093608 | 0.759370 | 0.446460 | Marti Neuron |
| 6 | 951094050 | 4.207487 | -1 | 0.529467 | 951092695 | 0.848769 | 0.435930 | Other Neurons |
| 7 | 951092398 | 3.808470 | 1 | 0.655551 | 951092410 | 0.378662 | 0.428909 | Other Neurons |
| 8 | 951093928 | 4.130337 | 1 | 0.750732 | 951098242 | 0.877427 | 0.418735 | Other Neurons |
| 9 | 951094081 | 3.507181 | 1 | 0.848041 | 951098242 | 0.603391 | 0.403576 | Other Neurons |
| 10 | 951093608 | 4.435348 | -1 | 0.736594 | 951093821 | 2.418180 | 0.400721 | Other Neurons |
| 11 | 951093725 | 2.987012 | 1 | 0.675827 | 951092940 | 0.228417 | 0.364932 | Other Neurons |
| 12 | 951093662 | 3.655037 | 1 | 0.602763 | 951094559 | 0.284722 | 0.353167 | Other Neurons |
| 13 | 951099306 | 3.979247 | 1 | 0.797598 | 951094081 | 0.391902 | 0.348804 | Other Neurons |
| 14 | 951092949 | 4.111789 | -1 | 0.611074 | 951093608 | 1.177489 | 0.323257 | Other Neurons |
| 15 | 951092940 | 4.349405 | -1 | 0.675827 | 951093725 | 2.236061 | 0.320312 | Other Neurons |
| 16 | 951093283 | 4.400400 | -1 | 0.414326 | 951092437 | 1.576702 | 0.313524 | Other Neurons |
| 17 | 951094516 | 4.148279 | 1 | 0.474654 | 951092949 | 0.410199 | 0.295350 | Other Neurons |
| 18 | 951092973 | 4.301265 | -1 | 0.559254 | 951097547 | 1.121742 | 0.295022 | Other Neurons |
| 19 | 951092410 | 4.041892 | 1 | 0.655551 | 951092398 | 0.579999 | 0.293418 | Other Neurons |
| 20 | 951093349 | 3.616517 | 1 | 0.482599 | 951092410 | 0.193406 | 0.291107 | Other Neurons |
| 21 | 951098928 | 4.362596 | -1 | 0.503125 | 951092437 | 1.434995 | 0.286298 | Other Neurons |
| 22 | 951092924 | 3.903790 | 1 | 0.252571 | 951098384 | 0.158347 | 0.279925 | Other Neurons |
| 23 | 951093472 | 4.131069 | 1 | 0.497556 | 951093662 | 0.828043 | 0.279441 | Other Neurons |
| 24 | 951092961 | 4.271673 | -1 | 0.478252 | 951092983 | 0.805439 | 0.276420 | Other Neurons |
| 25 | 951093688 | 4.336214 | 1 | 0.672685 | 951093608 | 0.824899 | 0.271322 | Other Neurons |
| 26 | 951097826 | 4.057939 | 1 | 0.703051 | 951098560 | 0.241768 | 0.260815 | Other Neurons |
| 27 | 951098332 | 3.943135 | 1 | 0.363750 | 951097653 | 0.279123 | 0.259099 | Other Neurons |
| 29 | 951097708 | 3.697292 | 1 | 0.559016 | 951093846 | 0.174054 | 0.239430 | Other Neurons |
| 30 | 951092450 | 3.951929 | 1 | 0.677264 | 951093846 | 0.244408 | 0.235229 | Other Neurons |
| 31 | 951092983 | 4.170997 | 1 | 0.539202 | 951097919 | 0.472773 | 0.228897 | Other Neurons |
| 32 | 951093484 | 3.909350 | 1 | 0.397126 | 951094529 | 0.170366 | 0.224352 | Other Neurons |
| 33 | 951097504 | 4.165286 | 1 | 0.565490 | 951092450 | 0.556476 | 0.223981 | Other Neurons |
| 34 | 951097919 | 4.069588 | 1 | 0.539202 | 951092983 | 0.254956 | 0.217594 | Other Neurons |
| 35 | 951097653 | 3.692260 | 1 | 0.617355 | 951098360 | 0.231239 | 0.213991 | Other Neurons |
| 36 | 951097934 | 3.846258 | 1 | 0.691230 | 951093125 | 0.196590 | 0.204507 | Other Neurons |
| 37 | 951099320 | 3.872860 | 1 | 0.462756 | 951098242 | 0.191630 | 0.204198 | Other Neurons |
| 38 | 951098105 | 4.005780 | 1 | 0.485009 | 951094081 | 0.360385 | 0.202622 | Other Neurons |
| 39 | 951098075 | 3.906645 | 1 | 0.380846 | 951093928 | 0.310396 | 0.201916 | Other Neurons |
| 40 | 951097857 | 3.859465 | 1 | 0.582236 | 951093662 | 0.156676 | 0.193375 | Other Neurons |
| 41 | 951098773 | 4.007322 | 1 | 0.590119 | 951093125 | 0.371217 | 0.192125 | Other Neurons |
| 42 | 951092303 | 4.087908 | -1 | 0.416028 | 951093160 | 0.552178 | 0.186200 | Other Neurons |
| 43 | 951092437 | 4.423117 | -1 | 0.515452 | 951097547 | 1.912515 | 0.181156 | Other Neurons |
| 45 | 951097869 | 3.996161 | -1 | 0.308719 | 951092594 | 0.522796 | 0.175326 | Other Neurons |
| 46 | 951093069 | 3.957978 | 1 | 0.428547 | 951093125 | 0.173575 | 0.174352 | Other Neurons |
| 47 | 951098265 | 3.896537 | 1 | 0.365088 | 951098773 | 0.125817 | 0.170916 | Other Neurons |
| 49 | 951099491 | 3.959182 | 1 | 0.403545 | 951092303 | 0.179624 | 0.161948 | Other Neurons |
| 51 | 951093160 | 3.919836 | 1 | 0.416028 | 951092303 | 0.127743 | 0.140241 | Other Neurons |
| 52 | 951098560 | 3.918091 | 1 | 0.703051 | 951097826 | 0.307136 | 0.127194 | Other Neurons |
top_oracle_neurons.to_csv('top_oracle_neurons.csv', index=False)
In order to select the most informative neurons, we developed an oracle correlation method which analyses a neuron's response across different images. It compares a random response to the average response on other trials for the same image and penalises neurons with limited image data. \
To ensure reliability, we developed a second function which runs this process 50 times (i.e., mean number of trials) and calculates average scores. \ We focused on high-performing neurons with strong average oracle correlations (above 60th percentile) and high variance (above 50th percentile).
In the following we define some auxiliary functions necessary for later computations of our Convolutional Neural Network. A more detailed description of the model can be found below.
def laplace():
return np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]]).astype(np.float32)[None, None, ...]
def positive(weight, cache=None):
weight.data *= weight.data.ge(0).float()
return cache
class Pyramid(nn.Module):
_filter_dict = {
'gauss5x5': np.float32([
[0.003765, 0.015019, 0.023792, 0.015019, 0.003765],
[0.015019, 0.059912, 0.094907, 0.059912, 0.015019],
[0.023792, 0.094907, 0.150342, 0.094907, 0.023792],
[0.015019, 0.059912, 0.094907, 0.059912, 0.015019],
[0.003765, 0.015019, 0.023792, 0.015019, 0.003765]]),
'gauss3x3': np.float32([
[1 / 16, 1 / 8, 1 / 16],
[1 / 8, 1 / 4, 1 / 8],
[1 / 16, 1 / 8, 1 / 16]]
),
'laplace5x5': np.outer(np.float32([1, 4, 6, 4, 1]), np.float32([1, 4, 6, 4, 1])) / 256,
}
def __init__(self, scale_n=4, type='gauss5x5', downsample=True, _skip_upsampling=False):
"""
Setup Laplace image pyramid
Args:
scale_n: number of Laplace pyramid layers to construct
type: type of Gaussian filter used in pyramid construction. Valid options are: 'gauss5x5', 'gauss3x3', and 'laplace5x5'
downsample: whether to downsample the image in each layer. Defaults to True
"""
super().__init__()
self.type = type
self.downsample = downsample
self._skip_upsampling = _skip_upsampling
h = self._filter_dict[type]
self.register_buffer('filter', torch.from_numpy(h))
self.scale_n = scale_n
self._kern = h.shape[0]
self._pad = self._kern // 2
def lap_split(self, img):
N, c, h, w = img.size()
filter = self.filter.expand(c, 1, self._kern, self._kern).contiguous()
# the necessary output padding depends on even/odd of the dimension
output_padding = (h + 1) % 2, (w + 1) % 2
smooth = F.conv2d(img, filter, padding=self._pad, groups=c)
if self.downsample:
lo = smooth[:, :, ::2, ::2]
if self._skip_upsampling:
# Technically incorrect implementation of the Laplace pyramid
lo2 = smooth
else:
lo2 = 4 * F.conv_transpose2d(lo, filter, stride=2, padding=self._pad, output_padding=output_padding,
groups=c)
else:
lo = lo2 = smooth
hi = img - lo2
return lo, hi
def forward(self, img):
levels = []
for i in range(self.scale_n):
img, hi = self.lap_split(img)
levels.append(hi)
levels.append(img)
return levels
def __repr__(self):
return "Pyramid(scale_n={scale_n}, padding={_pad}, downsample={downsample}, type={type})".format(
**self.__dict__)
class SpatialTransformerPyramid2d(nn.Module):
def __init__(self, in_shape, outdims, scale_n=4, positive=False, bias=True,
init_range=.1, downsample=True, _skip_upsampling=False, type=None):
super().__init__()
self.in_shape = in_shape
c, w, h = in_shape
#print(f'c: {c}')
self.outdims = outdims
self.positive = positive
self.gauss_pyramid = Pyramid(scale_n=scale_n, downsample=downsample, _skip_upsampling=_skip_upsampling, type=type)
self.grid = Parameter(torch.Tensor(1, outdims, 1, 2))
self.features = Parameter(torch.Tensor(1, c * (scale_n + 1), 1, outdims))
#print(f'self.features first: {self.features.size()}')
if bias:
bias = Parameter(torch.Tensor(outdims))
self.register_parameter('bias', bias)
else:
self.register_parameter('bias', None)
self.init_range = init_range
self.initialize()
def initialize(self):
self.grid.data.uniform_(-self.init_range, self.init_range)
self.features.data.fill_(1 / self.in_shape[0])
if self.bias is not None:
self.bias.data.fill_(0)
def group_sparsity(self, group_size):
f = self.features.size(1)
n = f // group_size
ret = 0
for chunk in range(0, f, group_size):
ret = ret + (self.features[:, chunk:chunk + group_size, ...].pow(2).mean(1) + 1e-12).sqrt().mean() / n
return ret
def feature_l1(self, average=True):
if average:
return self.features.abs().mean()
else:
return self.features.abs().sum()
def neuron_layer_power(self, x, neuron_id):
if self.positive:
positive(self.features)
self.grid.data = torch.clamp(self.grid.data, -1, 1)
N, c, w, h = x.size()
m = self.gauss_pyramid.scale_n + 1
feat = self.features.view(1, m * c, self.outdims)
print(f'self.features second: {feat.size()}')
y = torch.cat(self.gauss_pyramid(x), dim=1)
y = (y * feat[:, :, neuron_id, None, None]).sum(1)
return y.pow(2).mean()
def forward(self, x, shift=None):
if self.positive:
positive(self.features)
self.grid.data = torch.clamp(self.grid.data, -1, 1)
N, c, w, h = x.size()
m = self.gauss_pyramid.scale_n + 1
feat = self.features.view(1, m * c, self.outdims)
if shift is None:
grid = self.grid.expand(N, self.outdims, 1, 2)
else:
grid = self.grid.expand(N, self.outdims, 1, 2) + shift[:, None, None, :]
pools = [F.grid_sample(xx, grid, align_corners=True) for xx in self.gauss_pyramid(x)]
y = torch.cat(pools, dim=1).squeeze(-1)
y = (y * feat).sum(1).view(N, self.outdims)
if self.bias is not None:
y = y + self.bias
return y
def __repr__(self):
c, w, h = self.in_shape
r = self.__class__.__name__ + \
' (' + '{} x {} x {}'.format(c, w, h) + ' -> ' + str(self.outdims) + ')'
if self.bias is not None:
r += ' with bias'
for ch in self.children():
r += ' -> ' + ch.__repr__() + '\n'
return r
class Laplace(nn.Module):
"""
Laplace filter for a stack of data.
"""
def __init__(self, padding=0):
super().__init__()
self._padding = padding
self.register_buffer('filter', torch.from_numpy(laplace()))
def forward(self, x):
return F.conv2d(x, self.filter, padding=self._padding, bias=None)
class LaplaceL2(nn.Module):
"""
Laplace regularizer for a 2D convolutional layer.
"""
def __init__(self, padding=0):
super().__init__()
self.laplace = Laplace(padding=padding)
def forward(self, x, weights=None):
ic, oc, k1, k2 = x.size()
if weights is None:
weights = 1.0
return (self.laplace(x.view(ic * oc, 1, k1, k2)).view(ic, oc, k1, k2).pow(2) * weights).mean() / 2
class Core:
def initialize(self):
#log.info('Not initializing anything')
pass
def __repr__(self):
s = super().__repr__()
s += ' [{} regularizers: '.format(self.__class__.__name__)
ret = []
for attr in filter(lambda x: 'gamma' in x or 'skip' in x, dir(self)):
ret.append('{} = {}'.format(attr, getattr(self, attr)))
return s + '|'.join(ret) + ']\n'
class Core2d(Core):
def initialize(self, cuda=False):
self.apply(self.init_conv)
if cuda:
self = self.cuda()
@staticmethod
def init_conv(m):
if type(m) == nn.Conv2d:
nn.init.xavier_normal_(m.weight.data)
if m.bias is not None:
nn.init.constant_(m.bias.data, 0)
#def init_conv(m):
# if isinstance(m, nn.Conv2d):
# nn.init.xavier_normal_(m.weight.data)
# if m.bias is not None:
# m.bias.fill_(0)
# ---------------------- Identity Core ---------------------------
class IdentityCore(nn.Module, Core2d):
def __init__(self, **kwargs):
#log.info('Ignoring input {} when creating {}'.format(repr(kwargs), self.__class__.__name__))
super().__init__()
def forward(self, x):
return x
def regularizer(self):
return 0.0
# ---------------------- Conv2d Cores with Pooling -----------------------------
class Stacked2dCore_pooling(Core2d, nn.Module):
def __init__(self, input_channels, hidden_channels, input_kern, hidden_kern, pooling, hidden_pooling, layers=3,
gamma_hidden=0, gamma_input=0., skip=0, final_nonlinearity=True, bias=False, skip_nonlin=False,
momentum=0.1, pad_input=True, batch_norm=True, laplace_padding=0, laplace_weights_fn=None,**kwargs):
#log.info('Ignoring input {} when creating {}'.format(repr(kwargs), self._class.name_))
super().__init__()
if skip_nonlin:
print('Skip non-linearity')
self._input_weights_regularizer = LaplaceL2(padding=laplace_padding)
self.layers = layers
self.gamma_input = gamma_input
self.gamma_hidden = gamma_hidden
self.input_channels = input_channels
self.hidden_channels = hidden_channels
self.skip = skip
self.features = nn.Sequential()
self.pooling = pooling
self.hidden_pooling = hidden_pooling
# --- first layer
layer = OrderedDict()
layer['conv'] = \
nn.Conv2d(input_channels, hidden_channels, input_kern,
padding=input_kern // 2 if pad_input else 0, bias=bias)
layer['pool'] = nn.AvgPool2d(pooling)
if batch_norm:
layer['norm'] = nn.BatchNorm2d(hidden_channels, momentum=momentum)
if (not skip_nonlin) and (layers > 1 or final_nonlinearity):
layer['nonlin'] = nn.ELU(inplace=True)
self.features.add_module('layer0', nn.Sequential(layer))
# --- other layers
for l in range(1, self.layers):
layer = OrderedDict()
layer['conv'] = \
nn.Conv2d(hidden_channels if not skip > 1 else min(skip, l) * hidden_channels,
hidden_channels, hidden_kern,
padding=hidden_kern // 2, bias=bias)
layer['pool'] = nn.AvgPool2d(hidden_pooling)
if batch_norm:
layer['norm'] = nn.BatchNorm2d(hidden_channels, momentum=momentum)
if (not skip_nonlin) and (final_nonlinearity or l < self.layers - 1):
layer['nonlin'] = nn.ELU(inplace=True)
self.features.add_module('layer{}'.format(l), nn.Sequential(layer))
self.apply(self.init_conv)
if laplace_weights_fn is not None:
_, _, h, w = self.features[0].conv.weight.size()
dist_grid = torch.sqrt(torch.linspace(-1, 1, h)[:, None].pow(2) + torch.linspace(-1, 1, w).pow(2))
self.register_buffer('laplace_weights', laplace_weights_fn(dist_grid))
else:
self.laplace_weights = None
def forward(self, input_):
ret = []
for l, feat in enumerate(self.features):
do_skip = l >= 1 and self.skip > 1
input_ = feat(input_ if not do_skip else torch.cat(ret[-min(self.skip, l):], dim=1))
ret.append(input_)
return torch.cat(ret, dim=1)
def laplace(self):
return self._input_weights_regularizer(self.features[0].conv.weight, weights=self.laplace_weights)
def group_sparsity(self):
ret = 0
for l in range(1, self.layers):
ret = ret + self.features[l].conv.weight.pow(2).sum(3, keepdim=True).sum(2, keepdim=True).sqrt().mean()
return ret / ((self.layers - 1) if self.layers > 1 else 1)
def regularizer(self):
return self.group_sparsity() * self.gamma_hidden + self.gamma_input * self.laplace()
@property
def outchannels(self):
return len(self.features) * self.hidden_channels
#951092594
class CustomDataset(Dataset):
def __init__(self, root_dir, dataframe, transform=None, target_size=(60, 96)):
self.root_dir = root_dir
self.dataframe = dataframe
self.transform = transform
self.target_size = target_size
# Filter out non-image files and hidden files
self.image_names = [f for f in os.listdir(self.root_dir) if os.path.isfile(os.path.join(self.root_dir, f)) and not f.startswith('.')]
self.responses = self.dataframe['951092075'].tolist()
def __len__(self):
return len(self.image_names)
def __getitem__(self, idx):
image_name = self.image_names[idx]
image_path = os.path.join(self.root_dir, image_name)
# Load image
image = Image.open(image_path).convert('L') # Convert to grayscale
# Resize image
image = Resize(self.target_size)(image)
# Apply additional transformations
if self.transform:
image = self.transform(image)
# Get corresponding neuronal response
response = self.responses[idx]
return image, response
class CustomCNN_2_pooling(nn.Module):
def __init__(self, input_channels):
super(CustomCNN_2_pooling, self).__init__()
# Define the CNN core with three layers
self.conv_core = Stacked2dCore_pooling(input_channels= input_channels,
hidden_channels=64,
input_kern=4,
hidden_kern=3,
pooling = 2,
hidden_pooling=1,
layers=3,
gamma_hidden=0.1,
gamma_input=0.1,
skip=1,
final_nonlinearity=True,
bias=False,
skip_nonlin=False,
momentum=0.1,
pad_input=True,
batch_norm=True,
laplace_padding=0, stride = 3)
# Calculate the number of output channels from the CNN core
self.outchannels = self.conv_core.outchannels
#self.readout = CustomReadout(self.outchannels)
self.readout = SpatialTransformerPyramid2d(in_shape=(self.outchannels, 96, 60),
type = 'gauss5x5',
scale_n=3,
# outchannels = 256,
outdims=1,
positive=True)
def forward(self, x):
# Forward pass through the convolutional core
features = self.conv_core(x)
outputs = self.readout(features)
return outputs
The fundamental architecture comprises 3 blocks, each followed by a neuron-specialized readout layer. More specifically, the structure is as follows.
Convolutional Layer: Utilizes a kernel size tailored to match available Gabor filters, facilitating an effective comparison between stimuli.
Batch Normalization Layer: Essential for preventing overfitting by standardizing and normalizing the inputs of each layer.
Non-linearity Activation: Implemented using the Exponential Linear Unit (ELU) function, defined as:
This structured approach optimizes the network's performance by enabling effective feature extraction, regularization, and computational efficiency.
The readout layer synthesizes the features extracted by the CNN into a single scalar representing the average spike activity of a neuron. This layer is modeled as an affine transformation of the features, followed by a non-linearity. To ensure the output remains positive, the activation function employed is $ELU(x)+1$.
For further insights into the entire network structure, including the readout layer, refer to the following paper:
Walker, E.Y., Sinz, F.H., Cobos, E. et al. "Inception loops discover what excites neurons most using deep predictive models." Nature Neuroscience, vol. 22, pp. 2060–2065, 2019.
This approach encapsulates the essence of the network's functionality, providing a mechanism for interpreting neuron excitations through deep predictive models.
In this section, we proceed by generating the Most Excitatory Input for the selected neurons. We will exploit the convolutional neural network defined and trained in the previous section so that, by identifying which features of the stimuli (such as edges, colors, orientations, movements, etc.) trigger the highest firing rates, it is now possible can construct a stimulus that combines these features in a way that is most likely to activate the neuron maximally.
We decided to set up a new dataset by merging the training and the test set of the CNN. We now have a dataset that has the 118 natural grayscale images as indices and the neurons of the session as columns. The entries are torch tensors of reduced dimension (1, 60, 96) and encode the images' pixel intensity.
# Here I import a dataset containing all 118 natural scenes of a
mei_dir = '/content/drive/MyDrive/MEI testing/Images_final'
dataframe_mei = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/dataset_for_agne_mean.csv')
dataframe_mei.sort_values(by='frame', inplace=True)
dataframe_mei.set_index('frame', inplace=True, drop=True)
dataframe_mei = dataframe_mei.drop(dataframe_mei.index[0])
# Here I let the images be reduced in 60X96 format
mei_dataset = CustomDataset(mei_dir, dataframe_mei, transform=transforms.ToTensor())
Below, we provide an example of one of the images or our dataset. It is possible to notice that it is more blurred than the original, and that the edges are less sharp than a high-resolution picture.
a = np.random.randint(0,117)
plt.imshow(mei_dataset[a][0].squeeze().numpy(), cmap='gray')
plt.axis('off')
plt.show()
Then, we load the model we stored in the following checkpoint file and set it to evaluation mode. \ This corresponds to the model used to predict the responses for the neuron $\text{neuron_id} = '951097947'$, which is associated with the highest recorded accuracy value ($59$%).
neuron_id = '951097947'
# Define the path to the checkpoint file
model_path = '/content/drive/MyDrive/Neuroscience Project/checkpoints_final/checkpoint_epoch_951097947_size_gabor_kernel_59_accuracy.pth'
# Initialize the model
model = CustomCNN_2_pooling(1) # Replace YourModelClass with the class name of your model
# Load the checkpoint
checkpoint = torch.load(model_path)
# Load the model state dictionary
model.load_state_dict(checkpoint['model_state_dict'])
# Optional: If you have an optimizer saved in the checkpoint and want to load it
# optimizer = YourOptimizerClass(model.parameters(), lr=your_learning_rate) # Replace YourOptimizerClass with the class name of your optimizer
# optimizer.load_state_dict(checkpoint['optimizer_state_dict'])
# Access other saved information
epoch = checkpoint['epoch']
train_loss = checkpoint['train_loss']
val_loss = checkpoint['val_loss']
train_correlation = checkpoint['train_correlation']
val_correlation = checkpoint['val_correlation']
# Optionally, you can set the model to evaluation mode if needed
model.eval()
Given the extended dataset we constructed, we now compute the average contrast $C^*$ (defined as the mean standard deviation) of all images. This variable represents the target constant value of the MEI that will be generated, and will be used later on in the optimization to construct the penalty term of our loss function to optimize. \ The main purpose of $C^*$ is to help adjust the contrast of the MEI at each iteration, in order produce a stimulus that keeps roughly the same characteristics as the dataset the CNN was trained on.
# Initialization
c_star_std = 0
for i in range(len(mei_dataset)):
#std of images
std = torch.std(mei_dataset[i][0])
c_star_std += std
# Average over all images
c_star_std /= len(mei_dataset)
# Put in into number format
c_star_std = c_star_std.item()
print(f"Mean standard deviation for our dataset is {c_star_std}")
Mean standard deviation for our dataset is 0.2655758261680603
In order to proceed with the generation of the MEI, we follow the steps of the algorithm:
Repeat 2-5 for a preferred number of steps.\ Moreover, in order to keep the image's contrast aligned with respect to the average standard deviation of the natural scenes, we performed a L1-regularization in the loss function: $ \textit{L} = r - λ⋅|c^* - c| $
We proceed by defininig the functions that will be needed for the actual computation: \ $\text{grad}$: extracts the gradient from the $\text{backward()}$ pass of the CNN; \ $\text{fft_smooth}$: performs a Fast Fourier transform, which helps in suppressing high frequency noises making the variables smoother for optimization; \ $\text{scale}$: used here to scale gradients in the range $[0,1]$.
# Grad computation function
def grad(model, img, i, r, stds):
# Set the model to evaluation mode (no weights update)
model.eval()
# Contrast regularization parameter
param = 2
# Ensure image is the correct dimension
img = img.view(1, 1, 60, 96)
# Ensure tracking of gradients for the input image
img.requires_grad_(True)
# Forward pass through the model
outputs = model(img)
output = outputs.sum()
# Save statistics
std = torch.std(img)
r.append(output)
stds.append(std)
# Regularize with respect to the computed c_star_std
output -= param * torch.abs(torch.tensor(c_star_std) - std)
# Print to see how the optimization is working
if i % 10 == 0:
print(f"Loss to minimize is {output} in step {i} with std {std}")
# Clear previous gradients
model.zero_grad()
# Backpropagate to compute gradients w.r.t. input images
output.backward()
# Get the gradients of the image and reshape to [H, W]
image_gradients = img.grad.detach().cpu().squeeze()
return image_gradients, output #output [60,96] torch tensor
# Define fast fourier transform
def fft_smooth(grad, factor= 1/4):
if factor == 0:
return grad.numpy() # Directly return the numpy array if no smoothing is to be applied
# Ensure grad is a torch tensor
if isinstance(grad, np.ndarray):
grad = torch.from_numpy(grad).float()
# Obtain dimensions
h, w = grad.shape
# Adjust tw for the reduced dimensionality of the rFFT output
tw = torch.arange(0, w//2 + 1).float()
th = torch.arange(0, h).float()
# Compute the frequency weights
# Distance from origin in frequency space, avoiding division by zero
radius = torch.sqrt(th[:, None]**2 + tw[None, :]**2) + 1e-8
t = 1 / radius**factor
# Normalize the filter to preserve overall energy
F = t / t.mean()
# Compute the FFT of the gradient
grad_fft = torch.fft.rfftn(grad, dim=(0, 1))
# Apply the filter in the frequency domain
filtered_fft = grad_fft * F
# Compute the inverse FFT to transform back to the spatial domain
filtered_grad = torch.fft.irfftn(filtered_fft, s=(h, w), dim=(0, 1))
# Convert back to numpy for any subsequent processing
return filtered_grad.numpy() #[60,96] numpy array
# Scale in a preferred interval [c,d]
def scale(img, c = 0,d = 1):
# Quantities for scaling
a = img.min()
b = img.max()
# Scaling formula
sc_matrix = c + (img - a) * (d - c) / (b - a)
return np.array(sc_matrix) # [60,96] numpy array
# Generate white Gaussian noise image
# scaled between [0, 1]
height, width = 60,96
wn = scale(np.random.normal(0, 1, size = (height, width)), 0, 1)
#display
plt.imshow(wn, cmap='gray')
plt.axis('off')
plt.show()
# Here we transform the white noise into a tensor [1,1,60,96] readable by our model
scaled_data_tensor = torch.unsqueeze(torch.tensor(normal_data), dim=0).float()
img_tens = torch.unsqueeze(scaled_data_tensor, 0)
print(img_tens.size())
torch.Size([1, 1, 60, 96])
Function containing everything needed to compute MEIs
The hyperparameters we aim to fine-tune to achieve an optimal performance are:
def generate_mei(wn, steps, sigma, dataset= mei_dataset, alpha=0.1, c_star= c_star_std):
# Here we transform the whote noise into a tensor readable by our model
scaled_data_tensor = torch.unsqueeze(torch.tensor(wn), dim=0).float()
img_tens = torch.unsqueeze(scaled_data_tensor, 0)
# Some statistics:
# Some Statistics
r,stds = [],[]
for i in range(1, steps+1):
# Compute the gradients for the imput image
gradients, output = grad(model, img_tens, i, r, stds)
response = output.item()
# Transform the input image as a numpy object for smoother computation
img_tens = img_tens.squeeze().detach().numpy()
# Fourier Transform to smooth gradients
fft_grad = fft_smooth(gradients)
# Do Gradient Ascent step
img_tens += alpha*fft_grad
# Perform a Gaussian Blur and reduce sigma at each step
sigma = sigma + ((0.001 - sigma)*i)/steps
img_tens = torch.tensor(scipy.ndimage.gaussian_filter(img_tens,sigma,order = 0))
# Scale to 0-1
img_tens = scale(img_tens,0,1)
# Bring it to tensor in order to re-initialize the for loop
img_tens = torch.unsqueeze(torch.unsqueeze(torch.tensor(img_tens), dim=0).float(), 0)
final_mei = img_tens.squeeze().numpy()
#final_response = response
return final_mei, r, stds
Now we can call
mei_59, r_59, std_59 = generate_mei(neuron_id)
and visualize the MEI.
mei_59, r_59, std_59 = generate_mei(wn, steps=1000, sigma=7)
Loss to minimize is 4.375064849853516 in step 10 with std 0.2515137791633606 Loss to minimize is 4.308009624481201 in step 20 with std 0.269641637802124 Loss to minimize is 4.273122310638428 in step 30 with std 0.2788967788219452 Loss to minimize is 4.265032768249512 in step 40 with std 0.2836967408657074 Loss to minimize is 4.283144474029541 in step 50 with std 0.2857617139816284 Loss to minimize is 4.344533920288086 in step 60 with std 0.28649529814720154 Loss to minimize is 4.481871128082275 in step 70 with std 0.2868041396141052 Loss to minimize is 4.77593469619751 in step 80 with std 0.28716763854026794 Loss to minimize is 5.128570079803467 in step 90 with std 0.2877623736858368 Loss to minimize is 5.434421539306641 in step 100 with std 0.2883569896221161 Loss to minimize is 5.6400580406188965 in step 110 with std 0.28283339738845825 Loss to minimize is 5.604672908782959 in step 120 with std 0.26160919666290283 Loss to minimize is 5.5008931159973145 in step 130 with std 0.2428705245256424 Loss to minimize is 5.414088726043701 in step 140 with std 0.2221423089504242 Loss to minimize is 5.378298759460449 in step 150 with std 0.19699597358703613 Loss to minimize is 5.35123348236084 in step 160 with std 0.17513230443000793 Loss to minimize is 5.330989837646484 in step 170 with std 0.15610796213150024 Loss to minimize is 5.315517902374268 in step 180 with std 0.1395726352930069 Loss to minimize is 5.303035736083984 in step 190 with std 0.12524454295635223 Loss to minimize is 5.292255401611328 in step 200 with std 0.11288338899612427 Loss to minimize is 5.282402992248535 in step 210 with std 0.10227282345294952 Loss to minimize is 5.27306604385376 in step 220 with std 0.09321900457143784 Loss to minimize is 5.264244079589844 in step 230 with std 0.08554375916719437 Loss to minimize is 5.256048202514648 in step 240 with std 0.07908265292644501 Loss to minimize is 5.248531341552734 in step 250 with std 0.07368449121713638 Loss to minimize is 5.241790294647217 in step 260 with std 0.06920875608921051 Loss to minimize is 5.235894203186035 in step 270 with std 0.06552410125732422 Loss to minimize is 5.230829238891602 in step 280 with std 0.0625109001994133 Loss to minimize is 5.2265305519104 in step 290 with std 0.060061972588300705 Loss to minimize is 5.222908020019531 in step 300 with std 0.05808281898498535 Loss to minimize is 5.219875812530518 in step 310 with std 0.05649103969335556 Loss to minimize is 5.217343330383301 in step 320 with std 0.05521591380238533 Loss to minimize is 5.215231418609619 in step 330 with std 0.054197508841753006 Loss to minimize is 5.213469505310059 in step 340 with std 0.05338575690984726 Loss to minimize is 5.211986064910889 in step 350 with std 0.05273943394422531 Loss to minimize is 5.21073579788208 in step 360 with std 0.05222510173916817 Loss to minimize is 5.2096781730651855 in step 370 with std 0.0518156997859478 Loss to minimize is 5.20878267288208 in step 380 with std 0.051489558070898056 Loss to minimize is 5.208019256591797 in step 390 with std 0.05122938007116318 Loss to minimize is 5.207365036010742 in step 400 with std 0.051021452993154526 Loss to minimize is 5.206799507141113 in step 410 with std 0.05085493251681328 Loss to minimize is 5.206311225891113 in step 420 with std 0.050721194595098495 Loss to minimize is 5.205882549285889 in step 430 with std 0.05061344429850578 Loss to minimize is 5.205508708953857 in step 440 with std 0.05052636191248894 Loss to minimize is 5.205178260803223 in step 450 with std 0.05045566335320473 Loss to minimize is 5.204885959625244 in step 460 with std 0.05039804428815842 Loss to minimize is 5.204625129699707 in step 470 with std 0.050350870937108994 Loss to minimize is 5.20439338684082 in step 480 with std 0.05031207948923111 Loss to minimize is 5.204188823699951 in step 490 with std 0.05028005689382553 Loss to minimize is 5.204003810882568 in step 500 with std 0.050253450870513916 Loss to minimize is 5.203835964202881 in step 510 with std 0.05023123323917389 Loss to minimize is 5.203686714172363 in step 520 with std 0.05021261051297188 Loss to minimize is 5.203553199768066 in step 530 with std 0.05019690468907356 Loss to minimize is 5.203433513641357 in step 540 with std 0.05018359050154686 Loss to minimize is 5.203327178955078 in step 550 with std 0.05017225444316864 Loss to minimize is 5.203231334686279 in step 560 with std 0.050162557512521744 Loss to minimize is 5.203144550323486 in step 570 with std 0.05015421658754349 Loss to minimize is 5.203067779541016 in step 580 with std 0.05014701932668686 Loss to minimize is 5.202995777130127 in step 590 with std 0.05014077574014664 Loss to minimize is 5.202930927276611 in step 600 with std 0.05013534426689148 Loss to minimize is 5.202873706817627 in step 610 with std 0.05013059824705124 Loss to minimize is 5.202821731567383 in step 620 with std 0.050126466900110245 Loss to minimize is 5.2027740478515625 in step 630 with std 0.05012283846735954 Loss to minimize is 5.202732563018799 in step 640 with std 0.05011964961886406 Loss to minimize is 5.202693939208984 in step 650 with std 0.05011684447526932 Loss to minimize is 5.202657699584961 in step 660 with std 0.05011436715722084 Loss to minimize is 5.202628135681152 in step 670 with std 0.05011218786239624 Loss to minimize is 5.20259952545166 in step 680 with std 0.050110235810279846 Loss to minimize is 5.202574729919434 in step 690 with std 0.05010852962732315 Loss to minimize is 5.202552795410156 in step 700 with std 0.05010700225830078 Loss to minimize is 5.2025322914123535 in step 710 with std 0.05010564625263214 Loss to minimize is 5.202513694763184 in step 720 with std 0.050104446709156036 Loss to minimize is 5.20249605178833 in step 730 with std 0.05010337382555008 Loss to minimize is 5.202479839324951 in step 740 with std 0.050102412700653076 Loss to minimize is 5.2024664878845215 in step 750 with std 0.05010156333446503 Loss to minimize is 5.20245361328125 in step 760 with std 0.050100814551115036 Loss to minimize is 5.2024431228637695 in step 770 with std 0.05010014772415161 Loss to minimize is 5.202430248260498 in step 780 with std 0.05009952932596207 Loss to minimize is 5.202418804168701 in step 790 with std 0.050098977982997894 Loss to minimize is 5.202408313751221 in step 800 with std 0.050098493695259094 Loss to minimize is 5.202400207519531 in step 810 with std 0.050098076462745667 Loss to minimize is 5.202393531799316 in step 820 with std 0.05009772256016731 Loss to minimize is 5.202387809753418 in step 830 with std 0.05009738728404045 Loss to minimize is 5.202383995056152 in step 840 with std 0.05009710043668747 Loss to minimize is 5.202380657196045 in step 850 with std 0.05009683221578598 Loss to minimize is 5.202377796173096 in step 860 with std 0.050096623599529266 Loss to minimize is 5.2023749351501465 in step 870 with std 0.05009642615914345 Loss to minimize is 5.202373027801514 in step 880 with std 0.05009624361991882 Loss to minimize is 5.202371120452881 in step 890 with std 0.050096094608306885 Loss to minimize is 5.202369689941406 in step 900 with std 0.05009595304727554 Loss to minimize is 5.202367782592773 in step 910 with std 0.050095826387405396 Loss to minimize is 5.202367305755615 in step 920 with std 0.050095703452825546 Loss to minimize is 5.202366828918457 in step 930 with std 0.05009560286998749 Loss to minimize is 5.202363967895508 in step 940 with std 0.05009549483656883 Loss to minimize is 5.202361583709717 in step 950 with std 0.05009540915489197 Loss to minimize is 5.2023606300354 in step 960 with std 0.0500953309237957 Loss to minimize is 5.202358722686768 in step 970 with std 0.050095248967409134 Loss to minimize is 5.202357769012451 in step 980 with std 0.050095196813344955 Loss to minimize is 5.202356815338135 in step 990 with std 0.05009514093399048 Loss to minimize is 5.20235538482666 in step 1000 with std 0.0500950887799263
plt.imshow(mei_59, cmap='gray')
plt.axis('off')
plt.show()
As we analyse how the predicted neuronal response to the MEI and its contrast evolve in time, we observe two insightful opposite behaviors, suggesting that the optimization is affecting the image during the process of MEI generation: the neuronal response that is predicted by the CNN at each iteration on the refined stimulus is increasing, and the standard deviation of the resulting image is decreasing.
This behavior indicates that the MEI is becoming more effective at exciting the neuron. This increase is expected because the optimization process is specifically designed to adjust the image to maximize the neuron's response. As MEIs are designed based on the unique response characteristics of the targeted neuron, the MEI evolves in such a way that it aligns more closely with the neuron's receptive field or the features that the neuron is most responsive to. Over iterations, the image is refined to enhance the features that maximize activation, thereby increasing its effectiveness in stimulating the neuron. \
The decreasing standard deviation in the image during this process can be attributed to several factors:
Convergence to Specific Features: As the optimization focuses on enhancing specific features that maximally activate the neuron, the image becomes more focused or 'purified' towards the most stimulating features. This results in a reduction of the overall variability (standard deviation), since the image loses its randomness and becomes more uniform.
Reduction of Noise: In the initial stages, the image may contain a lot of random noise (since it starts from white noise). As the optimization progresses, it selectively enhances certain features over others, often reducing the noise or less relevant features, leading to a decrease in the overall standard deviation of pixel intensities.
Normalization and Scaling Effects: Depending on how the image is being processed during the optimization (e.g., scaling, normalization) the process might naturally reduce the spread of pixel intensity values.
# Create a 1x2 grid of subplots
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
# Plot neural response
axs[0].plot(np.array([t.item() for t in r_59]), label='Neural Response')
axs[0].set_title(' Predicted Neural Response')
axs[0].set_xlabel('Steps')
axs[0].set_ylabel('Response Value')
# Plot standard deviation
axs[1].plot(np.array([t.item() for t in std_59]), label='Standard Deviation', color='orange')
axs[1].set_title('Standard Deviation')
axs[1].set_xlabel('Steps')
axs[1].set_ylabel('Standard Deviation')
# Add a legend
#axs[0].legend()
#axs[1].legend()
# Show the plots
plt.tight_layout()
plt.savefig("Statistics.png")
plt.show()
print(f'The maximal neuronal response is {r_59[-1]}')
print(f'The standard deviation of the MEI is {std_59[-1]}')
The maximal neuronal response is 5.203412055969238 The standard deviation of the MEI is 0.05023481696844101
While the increase in the neuronal response is a desired outcome indicating successful optimization towards an MEI, the decrease in standard deviation should be considered carefully:
Loss of Detail: If the standard deviation decreases significantly, it could imply that the image is losing detail. This might not necessarily be negative but is something to be aware of, especially if one may want to use this result to estimate the full spectrum of stimuli that affect the neuron.
Overfitting to Neuron: There's a possibility that the MEI is overfitting to the neuron's characteristics in the current setup, especially if the decrease in standard deviation suggests a loss of complexity in the image.
In order to ensure a successful result, it might be needed to compute the contrast of a masked version of the image, where only the circular area of the MEI is considered.
As it was previously stated, Most Excitatory Inputs are designed based on the unique response characteristics of the targeted neurons, combining features that maximize their activation. \ One may wonder whether the MEI depends on the image which the optimization procedure is initialized on. Therefore, we repeat the process feeding a black image (all pixels have intensity zero) to the function.
bl = np.zeros((60,96))
plt.imshow(bl, cmap='gray')
plt.axis('off')
plt.show()
mei_check, r_check, std_check = generate_mei(bl, steps=1000, sigma=7)
Predicted regularized response in step 10 is: 3.4271063804626465 with std 0.30462244153022766 Predicted regularized response in step 20 is: 3.278251886367798 with std 0.3353641629219055 Predicted regularized response in step 30 is: 3.264526844024658 with std 0.3388451039791107 Predicted regularized response in step 40 is: 3.2978897094726562 with std 0.3380120098590851 Predicted regularized response in step 50 is: 3.3261239528656006 with std 0.33806490898132324 Predicted regularized response in step 60 is: 3.362445116043091 with std 0.3379369378089905 Predicted regularized response in step 70 is: 3.4305732250213623 with std 0.33763590455055237 Predicted regularized response in step 80 is: 3.5854854583740234 with std 0.3372311294078827 Predicted regularized response in step 90 is: 3.7948713302612305 with std 0.3368215560913086 Predicted regularized response in step 100 is: 4.000244140625 with std 0.3364415764808655 Predicted regularized response in step 110 is: 4.246574401855469 with std 0.3277827501296997 Predicted regularized response in step 120 is: 4.559811115264893 with std 0.31225138902664185 Predicted regularized response in step 130 is: 4.872156620025635 with std 0.2989574074745178 Predicted regularized response in step 140 is: 5.156697750091553 with std 0.28635427355766296 Predicted regularized response in step 150 is: 5.407121181488037 with std 0.27407658100128174 Predicted regularized response in step 160 is: 5.553841590881348 with std 0.2559845447540283 Predicted regularized response in step 170 is: 5.575096607208252 with std 0.23008406162261963 Predicted regularized response in step 180 is: 5.586808204650879 with std 0.20652490854263306 Predicted regularized response in step 190 is: 5.590517997741699 with std 0.1852765530347824 Predicted regularized response in step 200 is: 5.586966037750244 with std 0.16623751819133759 Predicted regularized response in step 210 is: 5.577611923217773 with std 0.14879246056079865 Predicted regularized response in step 220 is: 5.564088344573975 with std 0.1334252953529358 Predicted regularized response in step 230 is: 5.534425735473633 with std 0.11956727504730225 Predicted regularized response in step 240 is: 5.503382205963135 with std 0.10764080286026001 Predicted regularized response in step 250 is: 5.472306728363037 with std 0.0974496528506279 Predicted regularized response in step 260 is: 5.442211151123047 with std 0.0888018012046814 Predicted regularized response in step 270 is: 5.413810729980469 with std 0.08151818066835403 Predicted regularized response in step 280 is: 5.387557029724121 with std 0.07543240487575531 Predicted regularized response in step 290 is: 5.363697528839111 with std 0.07038965076208115 Predicted regularized response in step 300 is: 5.3423261642456055 with std 0.06624595075845718 Predicted regularized response in step 310 is: 5.323410511016846 with std 0.06286844611167908 Predicted regularized response in step 320 is: 5.306825637817383 with std 0.06013665348291397 Predicted regularized response in step 330 is: 5.292392253875732 with std 0.057943157851696014 Predicted regularized response in step 340 is: 5.279891490936279 with std 0.05619363486766815 Predicted regularized response in step 350 is: 5.26910400390625 with std 0.054806943982839584 Predicted regularized response in step 360 is: 5.259816646575928 with std 0.05371422320604324 Predicted regularized response in step 370 is: 5.251826763153076 with std 0.052857816219329834 Predicted regularized response in step 380 is: 5.244958877563477 with std 0.05218995735049248 Predicted regularized response in step 390 is: 5.23906135559082 with std 0.05167163163423538 Predicted regularized response in step 400 is: 5.2339959144592285 with std 0.051271095871925354 Predicted regularized response in step 410 is: 5.229645729064941 with std 0.05096295475959778 Predicted regularized response in step 420 is: 5.225910186767578 with std 0.05072701722383499 Predicted regularized response in step 430 is: 5.222698211669922 with std 0.05054721608757973 Predicted regularized response in step 440 is: 5.219935894012451 with std 0.05041096732020378 Predicted regularized response in step 450 is: 5.217560768127441 with std 0.05030841380357742 Predicted regularized response in step 460 is: 5.215515613555908 with std 0.05023181438446045 Predicted regularized response in step 470 is: 5.213753700256348 with std 0.050175171345472336 Predicted regularized response in step 480 is: 5.212235450744629 with std 0.05013379827141762 Predicted regularized response in step 490 is: 5.21092414855957 with std 0.0501040555536747 Predicted regularized response in step 500 is: 5.209792137145996 with std 0.05008315667510033 Predicted regularized response in step 510 is: 5.208813190460205 with std 0.05006889998912811 Predicted regularized response in step 520 is: 5.2079668045043945 with std 0.05005962401628494 Predicted regularized response in step 530 is: 5.207233905792236 with std 0.050054050981998444 Predicted regularized response in step 540 is: 5.206599235534668 with std 0.05005118250846863 Predicted regularized response in step 550 is: 5.206047058105469 with std 0.05005024001002312 Predicted regularized response in step 560 is: 5.205570697784424 with std 0.05005071312189102 Predicted regularized response in step 570 is: 5.205156326293945 with std 0.05005211755633354 Predicted regularized response in step 580 is: 5.204796314239502 with std 0.0500541515648365 Predicted regularized response in step 590 is: 5.204482078552246 with std 0.050056569278240204 Predicted regularized response in step 600 is: 5.204209804534912 with std 0.050059184432029724 Predicted regularized response in step 610 is: 5.203973293304443 with std 0.0500619076192379 Predicted regularized response in step 620 is: 5.2037672996521 with std 0.05006461217999458 Predicted regularized response in step 630 is: 5.203588485717773 with std 0.050067245960235596 Predicted regularized response in step 640 is: 5.203433036804199 with std 0.050069764256477356 Predicted regularized response in step 650 is: 5.203299045562744 with std 0.05007214471697807 Predicted regularized response in step 660 is: 5.203179359436035 with std 0.05007435008883476 Predicted regularized response in step 670 is: 5.203074932098389 with std 0.050076380372047424 Predicted regularized response in step 680 is: 5.202983856201172 with std 0.05007825419306755 Predicted regularized response in step 690 is: 5.202906131744385 with std 0.05007999762892723 Predicted regularized response in step 700 is: 5.202838897705078 with std 0.05008157715201378 Predicted regularized response in step 710 is: 5.202776908874512 with std 0.050082966685295105 Predicted regularized response in step 720 is: 5.202723026275635 with std 0.05008423700928688 Predicted regularized response in step 730 is: 5.202676773071289 with std 0.05008537694811821 Predicted regularized response in step 740 is: 5.20263671875 with std 0.05008641257882118 Predicted regularized response in step 750 is: 5.202601909637451 with std 0.0500873401761055 Predicted regularized response in step 760 is: 5.20257043838501 with std 0.05008815973997116 Predicted regularized response in step 770 is: 5.202544212341309 with std 0.05008890479803085 Predicted regularized response in step 780 is: 5.202520847320557 with std 0.05008956044912338 Predicted regularized response in step 790 is: 5.202498912811279 with std 0.050090137869119644 Predicted regularized response in step 800 is: 5.202481269836426 with std 0.05009064823389053 Predicted regularized response in step 810 is: 5.202466011047363 with std 0.05009111762046814 Predicted regularized response in step 820 is: 5.202452659606934 with std 0.050091538578271866 Predicted regularized response in step 830 is: 5.202441215515137 with std 0.05009189993143082 Predicted regularized response in step 840 is: 5.2024312019348145 with std 0.05009223520755768 Predicted regularized response in step 850 is: 5.202423572540283 with std 0.05009251832962036 Predicted regularized response in step 860 is: 5.202415943145752 with std 0.05009276419878006 Predicted regularized response in step 870 is: 5.202408790588379 with std 0.050092991441488266 Predicted regularized response in step 880 is: 5.2024030685424805 with std 0.05009319633245468 Predicted regularized response in step 890 is: 5.20239782333374 with std 0.050093356519937515 Predicted regularized response in step 900 is: 5.202394008636475 with std 0.05009353160858154 Predicted regularized response in step 910 is: 5.202390670776367 with std 0.05009365826845169 Predicted regularized response in step 920 is: 5.20238733291626 with std 0.050093770027160645 Predicted regularized response in step 930 is: 5.202384948730469 with std 0.050093889236450195 Predicted regularized response in step 940 is: 5.202383041381836 with std 0.05009398236870766 Predicted regularized response in step 950 is: 5.202381134033203 with std 0.05009406432509422 Predicted regularized response in step 960 is: 5.202378273010254 with std 0.050094123929739 Predicted regularized response in step 970 is: 5.202377796173096 with std 0.05009419843554497 Predicted regularized response in step 980 is: 5.202375888824463 with std 0.05009424313902855 Predicted regularized response in step 990 is: 5.2023749351501465 with std 0.05009428411722183 Predicted regularized response in step 1000 is: 5.20237398147583 with std 0.050094325095415115
plt.imshow(mei_check, cmap='gray')
plt.axis('off')
plt.show()
In essence, MEIs are not random or natural stimuli but are instead crafted or optimized through computational methods to include the most excitatory features recognizable by the neuron.
As was mentioned earlier, one of the primary goals of creating MEIs is to understand what specific visual or sensory characteristics a neuron or group of neurons is most sensitive to. The conclusions for such a query could be provided by in-vivo validation; however, this approach is not feasible with the experimental setup we are working with. \
Comparing Most Excitatory Inputs (MEIs) with Gabor filters can indeed serve as a useful proxy for in vivo validation of neural models, especially in the context of visual neuroscience. \
Gabor filters are mathematical models that mimic the response of neurons in the mammalian visual cortex to specific visual stimuli, particularly those involving spatial frequency and orientation selectivity. \
We start with a dataset that follows the same setup as the pd.DataFrames used in the previous sections: we take the gabor filters' IDs as the DataFrame indices, and the neurons' IDs as the columns. This table allows to visualize the range of the spike counts triggered by the gabor filters for our selected neuron.
firing_rate_gabors = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/dataset_for_agne_gabor.csv')
# make it more readable
firing_rate_gabors = firing_rate_gabors.set_index(firing_rate_gabors['stimulus_presentation_id']).drop('stimulus_presentation_id', axis=1)
firing_rate_gabors.head()
| 951092050 | 951092075 | 951092303 | 951092369 | 951092398 | 951092410 | 951092437 | 951092450 | 951092475 | 951092488 | ... | 951098773 | 951098807 | 951098850 | 951098871 | 951098928 | 951099272 | 951099306 | 951099320 | 951099475 | 951099491 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| stimulus_presentation_id | |||||||||||||||||||||
| 1 | 1 | 4 | 5 | 2 | 1 | 4 | 0 | 1 | 1 | 5 | ... | 8 | 1 | 1 | 3 | 4 | 0 | 2 | 0 | 0 | 0 |
| 2 | 1 | 2 | 1 | 2 | 0 | 4 | 0 | 0 | 2 | 1 | ... | 10 | 2 | 2 | 7 | 5 | 0 | 3 | 0 | 0 | 0 |
| 3 | 2 | 6 | 2 | 2 | 0 | 8 | 3 | 0 | 1 | 1 | ... | 5 | 1 | 0 | 6 | 5 | 0 | 2 | 1 | 0 | 0 |
| 4 | 1 | 3 | 2 | 3 | 1 | 1 | 2 | 1 | 0 | 7 | ... | 6 | 2 | 0 | 3 | 0 | 0 | 2 | 1 | 0 | 0 |
| 5 | 3 | 4 | 1 | 3 | 0 | 1 | 15 | 3 | 2 | 0 | ... | 6 | 4 | 1 | 5 | 3 | 0 | 0 | 0 | 0 | 0 |
5 rows × 135 columns
firing_rate_gabors[neuron_id].min(), firing_rate_gabors[neuron_id].max()
(0, 8)
Using the DataFrame $\text{gabors_meta}$, it is possible to access the metadata about the Gabor filters in the same session of the Allen SDK Dataset.
gabors_meta = pd.read_csv('/content/drive/MyDrive/Neuroscience Project/gabors_meta.unknown', sep = '\t')
gabors_meta = gabors_meta.set_index('stimulus_presentation_id') # more readable:
gabors_meta.head()
| stimulus_block | start_time | stop_time | stimulus_name | x_position | size | phase | spatial_frequency | orientation | temporal_frequency | contrast | y_position | duration | stimulus_condition_id | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| stimulus_presentation_id | ||||||||||||||
| 1 | 0.0 | 84.942787 | 85.176306 | gabors | 20.0 | [20.0, 20.0] | [3644.93333333, 3644.93333333] | 0.08 | 90.0 | 4.0 | 0.8 | -10.0 | 0.233519 | 1 |
| 2 | 0.0 | 85.176306 | 85.426505 | gabors | 30.0 | [20.0, 20.0] | [3644.93333333, 3644.93333333] | 0.08 | 0.0 | 4.0 | 0.8 | 40.0 | 0.250199 | 2 |
| 3 | 0.0 | 85.426505 | 85.676704 | gabors | 40.0 | [20.0, 20.0] | [3644.93333333, 3644.93333333] | 0.08 | 45.0 | 4.0 | 0.8 | -30.0 | 0.250199 | 3 |
| 4 | 0.0 | 85.676704 | 85.926904 | gabors | 0.0 | [20.0, 20.0] | [3644.93333333, 3644.93333333] | 0.08 | 90.0 | 4.0 | 0.8 | -40.0 | 0.250199 | 4 |
| 5 | 0.0 | 85.926904 | 86.177122 | gabors | -40.0 | [20.0, 20.0] | [3644.93333333, 3644.93333333] | 0.08 | 45.0 | 4.0 | 0.8 | 10.0 | 0.250218 | 5 |
First, we observe that all gabors are presented to the mouse for a duration of approximately $0.25$ seconds, and that all gabors have roughly the same features. \
gabors_meta['duration'].value_counts().index
Index([0.2502017433212131, 0.2502182566787496, 0.2502157566788128,
0.2502067433212005, 0.2502132566787622, 0.2501967433212257,
0.2502042433212637, 0.2502207566788001, 0.2502157566786991,
0.2502182566788065,
...
0.2502017433213837, 0.2502257566789012, 0.2502092433212084,
0.2502232566786802, 0.2502142433212952, 0.2501992304961362,
0.2501942433212605, 0.2501917433212668, 0.2502132566787765,
0.2501917433208973],
dtype='float64', name='duration', length=230)
print('Sizes:', set(gabors_meta['size']))
print('Phases:', set(gabors_meta['phase']))
print('Spatial Frequencies:', set(gabors_meta['spatial_frequency']))
print('Temporal Frequencies:', set(gabors_meta['temporal_frequency']))
print('Contrast:', set(gabors_meta['contrast']))
Sizes: {'[20.0, 20.0]'}
Phases: {'[3644.93333333, 3644.93333333]'}
Spatial Frequencies: {0.08}
Temporal Frequencies: {4.0}
Contrast: {0.8}
The only parameters that differs across gabors are the orientation and position of the filter. However, we opted to exclude the position parameter when generating the most responsive gabor, as it was deemed irrelevant for the purpose of comparison with MEIs. Given the process used to generate MEIs, all filters are centered in our images. Consequently, we maintained the Gabor filters centered and focused solely on orientation. We observed only three balanced classes of orientations: $\{0°, 45°, 90°\}$.
gabors_meta['orientation'].value_counts()
orientation 90.0 1215 0.0 1215 45.0 1215 Name: count, dtype: int64
We now proceed with the extraction of the gabor filter that maximizes the chosen neuron's activation. \ We notice that, in case there exists more than a unique gabor that maximally excites the neuron, it can be useful to merge the first $k$ gabors that trigger the highest spike count. This approach not only leverages the collective excitatory potential of these top-performing Gabors but also addresses the constraints imposed by a limited parameter space. \
Therefore, the decision to merge these Gabors is driven by the need to introduce greater variability among the generated Gabors.
We define the following functions that allow to reconstruct the most excitatory gabor by taking the average of the $k=3$ gabors that trigger the highest firing rate:
## CREATE GABORS
def create_gabor(sz, phase, wavelength, orientation, sigma,
dy, dx):
""" Create a gabor patch (sinusoidal + gaussian).
Arguments:
height (int): Height of the image in pixels.
width (int): Width of the image in pixels.
phase (float): Angle at which to start the sinusoid in degrees.
wavelength (float): Wavelength of the sinusoid (1 / spatial frequency) in pixels.
orientation (float): Counterclockwise rotation to apply (0 is horizontal) in
degrees.
sigma (float): Sigma of the gaussian mask used in pixels
dy (float): Amount of translation in y (positive moves down) in pixels/height.
dx (float): Amount of translation in x (positive moves right) in pixels/height.
Returns:
Array of height x width shape with the required gabor.
"""
# Compute image size to avoid translation or rotation producing black spaces
padding = max(height, width)
imheight = height + 2 * padding
imwidth = width + 2 * padding
# we could have diff pad sizes per dimension = max(dim_size, sqrt((h/2)^2 + (w/2)^2))
# but this simplifies the code for just a bit of inefficiency
# Create sinusoid with right wavelength and phase
start_sample = phase
step_size = 360 / wavelength
samples = start_sample + step_size * np.arange(imheight)
samples = np.mod(samples, 360) # in degrees
rad_samples = samples * (np.pi / 180) # radians
sin = np.sin(rad_samples)
# Create Gabor by stacking the sinusoid along the cols
gabor = np.tile(sin, (imwidth, 1)).T
# Rotate around center
gabor = ndimage.rotate(gabor, orientation, reshape=False)
# Apply gaussian mask, its is normalized to take values in [0,1]
gaussy = signal.windows.gaussian(imheight, std=sigma)
gaussx = signal.windows.gaussian(imwidth, std=sigma)
mask = np.outer(gaussy, gaussx)
gabor = gabor * mask
# Translate (this is only approximate but it should be good enough)
if abs(dx) > 1 or abs(dy) > 1:
raise ValueError('Please express translations as factors of the height/width,'
'i.e, a number in interval [-1, 1] ')
dy = int(dy * height) # int is the approximation
dx = int(dx * width)
gabor = gabor[padding - dy: -padding - dy, padding - dx: -padding - dx]
if gabor.shape != (height, width):
raise ValueError('Dimensions of gabor do not match desired dimensions.')
return gabor.astype(np.float32)
def most_exc_gabors(neuron_id, df = firing_rate_gabors):
neuron_spks = df[neuron_id]
# Sort the indices based on the firing rates in descending order
sorted_indices = sorted(neuron_spks.index, key=lambda i: neuron_spks[i], reverse=True)
# Select the top 3 indices with the highest firing rates
top_3_indices = sorted_indices[:3]
return top_3_indices
## Computes the mean gabor as a numpy array
def mean_gabor(neuron_id, gab_meta = gabors_meta):
gabors_okay = gabors_meta.T.iloc[:,most_exc_gabors(neuron_id)].T
gabors_okay_ids = list(gabors_okay.index)
gabors = []
for i in gabors_okay_ids:
gabors.append(create_gabor(sz=(60,96), phase=3644.93333333, wavelength=np.pi, orientation=gabors_okay.loc[i,'orientation'], sigma=5,
dy=0, dx=0))
gabors_combined = np.stack(gabors, axis=0)
# Compute the mean across the arrays along axis 0
mean_gabors = np.mean(gabors_combined, axis=0)
# Scale in range [0, 1]
mean_gabors = scale(mean_gabors, 0, 1)
return mean_gabors
We can now visualize the gabor:
gabor_59 = mean_gabor(neuron_id)
plt.imshow(gabor_59, cmap='gray')
plt.axis('off')
plt.show()
Similarities in response: \ If the MEIs generated through computational models resemble the patterns produced by Gabor filters, this similarity can validate their biological plausibility. \
In the primary visual cortex, particularly in areas like V1, many neurons have receptive fields that respond optimally to specific orientations and spatial frequencies. These characteristics make Gabor patches ideal for studying and activating these neurons. However it is important to note that our study is limited to the fact that we cannot validate our findings through an in-vivo analysis, implying that MEIs and Gabor's response cannot be directly compared. Relevant research indicates that MEIs drive a higher response than Gabor patches for neurons in the V1 cortex. \
Where differences arise, it may indicate more complex or nonlinear processing than what Gabor filters can account for. Such differences could suggest additional layers of complexity in neural responses, especially in neurons with higher-order or more specialized functions. \
Validation of Neural Coding: By demonstrating that MEIs activate neurons in ways similar to Gabor filters, it is possible to infer that the computational model correctly identifies the types of stimuli to which real neurons are most responsive.
f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_59, cmap='gray')
axarr[0].axis('off')
axarr[1].imshow(gabor_59, cmap = 'gray');
axarr[1].axis('off')
plt.show()
When comparing Most Excitatory Inputs generated by deep learning models to traditional Gabor filters, it is not uncommon to find differences in appearance. While Gabor filters are specifically designed to model the spatial sensitivities of V1 simple cells, MEIs are derived from complex optimization processes that may incorporate a broader set of neuronal behaviors and response characteristics, including but not limited to orientation, spatial frequency, and phase. \
The images exhibit a subtle diagonal orientation; however, neither distinctly specifies a particular direction. This slight diagonal alignment suggests that the neuron is potentially responsive to stimuli oriented along this general diagonal axis, yet the precise optimal direction of activation remains undefined in both images.
Typically, Gabor filters are characterized by sinusoidal waves multiplied by a Gaussian envelope. This creates clear, periodic patterns across the filter and sharper edges than the MEI. Instead, in the MEI we notice the presence lower frequencies and wider patches of light and dark pixels.
These observations are insufficient to confirm the similarity between the two images or to fully elucidate the stimulus features to which the neuron is most selective. Consequently, we plan to explore some key metricsa deeper investigation is warranted to better understand the neuron's responsiveness. To this end, we plan to explore several key metrics, including luminance, contrast, 1-fold and 2-fold symmetry, and spatial frequency content.
Initially, it is necessary to create a circular mask that isolates only the central portion of the image where the Most Excitatory Input (MEI) is located. This region typically exhibits higher frequency features and a non-uniform background, distinguishing it from the peripheral areas of the image. This step is crucial because, by the end of the optimization process, it was observed that the overall contrast of the image was extremely low, likely due to the large uniform sections that did not contribute to neuron activation. To ensure that the computed values for luminance and contrast are not distorted by these uniform areas, they will be calculated specifically for the masked version of the image.
To estimate the radius of the circular mask, we simply look at the changes in pixel intensity along the horizontal axis at the level of the diameter.
# Let us check which is the most suitable radius for both gabors and MEIs
x = np.zeros(96)
for i in range(len(x)):
x[i] = mean_gabors[29][i]
plt.figure(figsize = (20,5))
plt.axvline(x=35, color='red', ls='--')
plt.axvline(x=61, color='red', ls='--')
plt.xticks([i for i in range(1,97)if i%2==0])
plt.plot(x)
[<matplotlib.lines.Line2D at 0x7edd589b3c40>]
MEI: we notice it is a bit shifted to the right (approx 5 pixels)
x = np.zeros(96)
for i in range(len(x)):
x[i] = final_mei[29][i]
plt.figure(figsize = (20,5))
plt.axvline(x=40, color='red', ls='--')
plt.axvline(x=66, color='red', ls='--')
plt.xticks([i for i in range(1,97)if i%2==0])
plt.plot(x)
[<matplotlib.lines.Line2D at 0x7edd588d2080>]
We observe a symmetric behavior and, in order to consider the variability in pixel values we find on the edges of the circle, we extend the radius to 12.
def create_circular_mask(h, w, center=None, radius=12):
if center is None: # use the middle of the image
center = (int(w/2), int(h/2))
if radius is None: # use the smallest distance between the center and image walls
radius = min(center[0], center[1], w-center[0], h-center[1])
Y, X = np.ogrid[:h, :w]
dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)
mask = dist_from_center <= radius
return mask
masked_59, lum_59 = luminance(mei_59, radius=12, center=(int(mei_59.shape[1]/2+3), int(mei_59.shape[0]/2)))
masked_gab_59, lum_gab_59 = luminance(gabor_59, radius=12, center=None)
f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(masked_59, cmap='gray')
axarr[1].imshow(masked_gab_59, cmap = 'gray');
The luminance value was calculated as the average pixel intensity within the spatial mask:
import math
def luminance(img, radius, center):
# apply mask:
h, w = img.shape[:2]
themask = create_circular_mask(h,w, radius=radius, center=center)
masked_img = img.copy()
masked_img[~themask] = 0.0
# compute luminance on masked img
lum=0
for x in masked_img:
lum += x.sum()
lum = lum / (math.pi*(radius**2))
return masked_img, lum
print(f'MEI luminance: {lum_59}')
print(f'Gabor luminance: {lum_gab_59}')
MEI luminance: 0.44963369290239985 Gabor luminance: 0.4097204913589961
From an initial visual inspection, the MEI appeared to exhibit higher luminance compared to the Gabor filter. However, when we delve into the metrics for a closer comparison, we find remarkable similarities between them. Notably, we observe that the luminance level that effectively activates the neuron is approximately $0.4$. This consistency across different visual tools underscores the specific luminance sensitivity of the neuron.
Contrast is a crucial factor in how neurons process visual stimuli, significantly more so than merely the lightness or darkness of an image. We thus proceed with the same localized approach as before:
print(f'C star: {c_star_std}')
C star: 0.2675706446170807
def masked_contrast(image, mask):
masked_values = image[mask] #Apply mask to the image
std_dev = np.std(masked_values) #standard deviation
return std_dev
print("Standard deviation of non-masked MEI:", mei_59.std())
print("Standard deviation of non-masked Gabor:", gabor_59.std())
Standard deviation of non-masked MEI: 0.050091017 Standard deviation of non-masked Gabor: 0.03668658
h, w = 60, 96
radius = 12
mask_mei_c = create_circular_mask(h, w, radius=radius, center= (int(w/2+3), int(h/2)))
masked_contrast_mei = masked_contrast(mei_59, mask_mei_c)
print("Masked Contrast MEI:", masked_contrast_mei)
mask_gab_c = create_circular_mask(h, w, radius=radius)
masked_contrast_gab = masked_contrast(gabor_59, mask_gab_c)
print("Masked Contrast Gabor:", masked_contrast_gab)
Masked Contrast MEI: 0.15114456 Masked Contrast Gabor: 0.13230467
We observe a significant difference betweeen the standard deviation of the unmasked image; moreover, we notice that the contrast values of masked images are very similar, and are not too far from the target contrast value $C^*$ identified as the mean contrast of our dataset.
Approach: \ To calculate folio (1-fold) and quarto (2-fold) , we first performed FFT and averaged the power spectrum axially with bin size = $\frac{\pi}{16}$ radians. To calculate the axial average, we converted pixel Cartesian coordinates to polar coordinates, defining the center of the image as the origin (0,0). Pixels were then binned according to their angle, rounded down. Next, we defined an n-fold symmetry index (SI) as $SI_n = \frac{|\sum_{\theta}P_{\theta}e^{i2n\theta}|}{\sum_\theta P_\theta}$ where $P_\theta$ is the average power at angle $\theta$ and the scalar multiplier 2 is due to the inherent point symmetry of the FFT power spectrum. SI is defined on the unit interval with 0 denoting a lack of symmetry and 1 being a fully n-fold symmetric image.
# GABOR
fft_gab = fftshift(fft2(gabor_59))
pow_spec_gab = (np.abs(fft_gab))**2
freq_x = np.fft.fftshift(np.fft.fftfreq(gabor_59.shape[1])) * gabor_59.shape[1]
freq_y = np.fft.fftshift(np.fft.fftfreq(gabor_59.shape[0])) * gabor_59.shape[0]
# MEI
fft_mei = fftshift(fft2(mei_59))
pow_spec_mei = (np.abs(fft_mei))**2
freq_x_mei = np.fft.fftshift(np.fft.fftfreq(mei_59.shape[1])) * mei_59.shape[1]
freq_y_mei = np.fft.fftshift(np.fft.fftfreq(mei_59.shape[0])) * mei_59.shape[0]
## function to compute axial average:
def axial_average(power, bin_size=1/16*np.pi):
height, width = power.shape
center = (int(width/2), int(height/2))
Y, X = np.indices((height, width))
#polar coordinates
R = np.sqrt((X - center[0])**2 + (Y - center[1])**2)
Theta = np.arctan2((Y - center[1]), (X - center[0]))
theta_bins = np.floor(Theta / bin_size).astype(int)
theta_bins += -theta_bins.min() # Shift scale to be positive
# Average by bins
radial_sum = np.bincount(theta_bins.ravel(), weights=power.ravel())
radial_count = np.bincount(theta_bins.ravel())
radial_mean = radial_sum / radial_count
return radial_mean, Theta, theta_bins
bin_size = 1/16*np.pi
P_theta, Theta, theta_bins = axial_average(pow_spec_gab)
P_theta_mei, Theta_mei, theta_bins_mei = axial_average(pow_spec_mei)
Here the 1-fold and 2-fold symmetry index follows:
def sym_index(P_theta, n):
angles = np.arange(len(P_theta)) * bin_size
ei2ntheta = np.exp(1j * 2 * n * angles)
numerator = np.abs(np.sum(P_theta * ei2ntheta))
denominator = np.sum(P_theta)
return numerator / denominator
# Calculate for 1-fold and 2-fold symmetry index
## GABOR
SI_1 = sym_index(P_theta, 1)
SI_2 = sym_index(P_theta, 2)
print("Gabor:")
print(f"1-Fold Symmetry Index: {SI_1}")
print(f"2-Fold Symmetry Index: {SI_2}\n")
## MEI
SI_1_mei = sym_index(P_theta_mei, 1)
SI_2_mei = sym_index(P_theta_mei, 2)
print("MEI:")
print(f"1-Fold Symmetry Index: {SI_1_mei}")
print(f"2-Fold Symmetry Index: {SI_2_mei}")
Gabor: 1-Fold Symmetry Index: 0.9992552829452915 2-Fold Symmetry Index: 0.9975716528845747 MEI: 1-Fold Symmetry Index: 0.9696106836897902 2-Fold Symmetry Index: 0.9845564551083268
The power spectrum helps in visualizing the energy distribution across different frequencies without considering the phase information. This simplification can be particularly useful for analyzing the overall impact of various frequency components, such as distinguishing between high and low-frequency trends in the data. \
In image processing, the power spectrum is useful to reveal details about the texture and orientation of objects within the image, as certain concentrations of energy in specific frequency bands can indicate the presence of particular textures or alignments.
f, axarr = plt.subplots(1,2, figsize=(30,10))
im0 = axarr[0].imshow(np.log1p(pow_spec_gab), extent=[freq_x.min(), freq_x.max(), freq_y.min(), freq_y.max()], cmap='gnuplot')
f.colorbar(im0, ax=axarr[0])
axarr[0].set_title("Log Power Spectrum of Gabor")
axarr[0].set_xlabel('horizontal freq')
axarr[0].set_ylabel('vertical freq')
im1 = axarr[1].imshow(np.log1p(pow_spec_mei), extent=[freq_x_mei.min(), freq_x_mei.max(), freq_y_mei.min(), freq_y_mei.max()], cmap='gnuplot');
f.colorbar(im1, ax=axarr[1])
axarr[1].set_title("Log Power Spectrum of MEI")
axarr[1].set_xlabel('horizontal freq')
axarr[1].set_ylabel('vertical freq')
plt.show()
We can now make a few observations:
Firstly, the zero frequency point at the center of the plot, often referred to as the DC (Direct Current) component in the context of Fourier analysis, represents the average or mean value of the entire image as it reflects the overall brightness or intensity. We observe that, as it was expected, in both the MEI and the Gabor, the zero frequency component is significant, due to a large uniform area in the background that does not vary much in intensity. \
While the Gabor's log power spectrum exhibits distinct peaks at high specific frequencies, the MEI's one displays a smooth, continuous distribution across the frequency range. \
This behavior is consistent with the structural difference between the two generated images: Gabor filters are designed to respond maximally at their characteristic frequencies, typically associated with the specific orientation and scale of the sinusoidal wave, which corresponds to distinct edges in the visual stimuli. For this reason, it is meaningful to notice that they display a very focused frequency response, showing a peak response at their characteristic frequency. \
In contrast, the dispersed frequency representation in the MEI's log power spectrum, with a significant contribution to the mid and low-frequency ranges, might suggest that the neuron is responsive to a broader variety of spatial frequencies, which aligns with processing less pronounced or smoother transitions in visual textures. \
Such pattern could also be justified by the slight presence of noise, which generally contributes to a more random spread across frequencies. \
Overall, since our MEI shows a more dispersed frequency representation, it could indicate that the neuron we examined responds not only to sharp, well-defined edges but also to edges with less contrast or those that are more diffusely defined. This can be a sign that the neuron processes information from a broader range of spatial frequencies.
Our analysis investigated the comparative efficacy of Most Excitatory Inputs (MEIs) and Gabor filters in maximizing neuronal responses for specific neurons. While MEIs are generally acknowledged as superior to Gabors, a direct comparison remains crucial for validating MEIs' specific advantages. Additionally, we want to investigate whether MEIs, optimized for visual stimuli response, exhibit parallels with Gabors, designed to evoke maximal responses in V1 neurons of the visual cortex. Our inquiry aimed to uncover shared characteristics between the two. Our findings reveal similarities between MEIs and Gabor patches, affirming the presence of common traits in driving neuronal response. Notably, research by Walker et al. (2019) suggests that MEIs elicit a heightened response compared to Gabor patches among V1 cortex neurons. However, it's essential to underscore that such analysis necessitates in vivo validation.
The limitations of our study stem from several factors:
Small Dataset: Both the number of recorded neurons and the quantity of natural images used for model training were limited. These limitations likely impacted the generalization capabilities of our model.
Nature of Deep Neural Networks: Deep neural networks often struggle to generalize beyond the typical statistics of their training data. Thus, it remains uncertain whether in silico synthesized MEIs can accurately predict in vivo responses, necessitating experimental verification of predictions derived from these models.
Computational Resources and Time Constraints: Resource limitations restricted the depth and breadth of our analysis, affecting both the complexity of our models and the scale of our simulations.
Image Dimensions and Quality: Due to limited computational resources, we resized the input images to expedite training, albeit at the expense of image quality. This compromise in image resolution hindered the model's ability to discern sharp edges and prominent patterns, thus impacting the learned characteristics of our model.
Limited Gabor Details: The parameter space for Gabor filters was constrained, primarily distinguishing only directions. With access to more diverse Gabors incorporating variations in phase, size, and eccentricity, we could have achieved a more precise comparison with our MEIs, thereby validating our analysis with greater accuracy. Moving forward, further analysis could involve:
Further Analysis
Expansion Beyond VISam Neurons: We focused exclusively on VISam neurons due to their prevalence in our available dataset, overlooking other areas of the visual cortex. Broadening our study to encompass neurons from diverse cortical regions could offer a more comprehensive insight into neuronal response patterns.
Exploration of Different Brain Regions: While detailed insights may not be expected, examining brain regions beyond V1 could offer additional perspectives on neural response predictability.
Comparison with Alternative Stimuli: Evaluating our model's performance against different types of stimuli could reveal its robustness and versatility.
Performance Comparison with Linear Models: Contrasting the performance of our deep neural network with that of linear models could elucidate the advantages and limitations of each approach.
Extension to Multi-Neuron Prediction: Expanding our CNN model to predict neuronal responses across multiple neurons simultaneously could enhance the sophistication and applicability of our simulations.
In the following section, we will only show other examples of Most Excitatory Inputs and Gabor filters corresponding to other neurons we selected in the first section of our analysis. \ In order to avoid repetitions, we decided to omit a detailed analysis for these neurons, as they are associated with a lower recorded accuracy value of the model used to predict their neuronal responses.
neuron_id = '951097919'
f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_46, cmap='gray')
axarr[1].imshow(gabor_46, cmap = 'gray');
neuron_id = '951092949'
f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_34, cmap='gray')
axarr[1].imshow(gabor_34, cmap = 'gray');
neuron_id = '951098928'
f, axarr = plt.subplots(1,2, figsize=(30,10))
axarr[0].imshow(mei_32, cmap='gray')
axarr[1].imshow(gabor_32, cmap = 'gray');